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PREFACE 


The  work  described  in  this  report  was  performed  under  Contract 
F33615-75-C-5054,  Task  Number  723105,  Work  Unit  7231-05-18  during  the 
period  from  December  1974  to  March  1976. 
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INTRODUCTION 


i 


Laving  bone  is  continually  undergoing  processes  of  growth,  reinforce- 
ment and  resorption.  These  processes  are  termed  collectively  "remodeling.  " 
The  remodeling  processes  in  living  bone  are  the  mechanisms  by  which  the 
bone  adapts  its  histological  structure  to  long-term  changes  in  loading.  Pro- 
longed straining  of  a bone  tends  to  make  the  bone  stiffer,  stronger  and  more 
dense.  This  is  why  a person  in  the  later  stages  of  recovery  from  a broken 
leg  is  encouraged  to  walk  on  his  healing  limb.  Conversely,  a living  bone  not 
subjected  to  its  accustomed  strain  level  will,  in  time,  become  more  compliant, 
weaker  and  more  porous.  An  example  of  this  occurs  when  men  are  subjected 
to  prolonged  weightlessness  such  as  in  the  recent  Skylab  experiments.  Re- 
searchers have  noted  that  after  these  astronauts  were  subjected  to  months  of 
weightlessness,  the  mineral  content  of  their  bones  had  decreased  significantly. 
They  estimate  that  for  longer  periods  of  weightlessness,  such  as  might  occur 
on  a flight  to  Mars,  this  skeletal  degeneration  (as  well  as  concomitant  atrophy 
in  the  cardiovascular  system)  might  become  dangerous.  Other  questions  of 
interest  involve  the  effects  of  prolonged  inactivity  on  the  strength  of  bones, 
and  therefore  on  the  physical  fitness  of  personnel.  For  this  reason  it  is  use- 
ful to  construct  a mathematical  model  describing  the  stress  adaptation  prop- 
erties of  living  bone  in  order  to  predict  what  exercise  regimen  might  be  most 
effective  in  stopping,  or  at  least  controlling,  the  skeletal  degeneration  brought 
on  by  weightlessness  and/or  physical  inactivity. 

This  report  describes  a mathematical  model  for  the  stress  adaptation 
of  cancellous  bone.  The  model  is  a generalization  of  the  mathematical  theory 
for  the  stress-adaptation  of  cortical  bone  originally  developed  by  the  author 
in  1971  and  first  presented  publicly  at  the  ASME  Biomechanics  Symposium 
in  Atlanta,  Georgia  in  June,  1973  (Hegedus,  1973).  More  detailed  develop- 
ment of  the  cortical  bone  model  can  be  found  in  the  .author's  dissertation 
(Hegedus,  1976)  and  in  two  other  papers  that  the  author  has  written  in  collab- 
oration with  Prof.  S.  C.  Cowin  (Cowin  and  Hegedus,  1976,  and  Hegedus  and 
Cowin,  1976).  The  development  of  the  theory  of  adaptive  elasticity  is  in  this 
paper. 

In  Section  I,  we  will  review  the  basic  anatomical  and  physiological 
properties  of  bone  which  bring  about  remodeling.  In  Section  II,  we  develop 
a general  theory  of  adaptive  elasticity  that  is  capable  of  describing  the 
trabecular  adaptation  of  cancellous  bone.  In  Section  III,  the  general  govern- 
ing equations  developed  in  the  previous  section  are  specialized  to  the  case  of 
two-dimensional  adaptation  to  plane  stress.  We  also  describe  and  cite  re- 
sults from  a computer  program  that  analyzes  problems  in  two-dimensional 
adaptive  elasticity.  Finally,  in  Section  IV,  the  experimental  procedure  is 
described  by  which  the  material  coefficients,  discussed  in  Section  II,  can  be 
measured  experimentally. 
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SECTION  I 

THE  BIOLOGICAL  PROPERTIES  OF  LIVING  BONE 

I.  1 INTRODUCTION 

In  this  Section  we  will  review  the  physiological  and  anatomical  properties 
of  bone  that  are  of  interest  in  the  present  study.  A more  detailed  description 
of  the  anatomy  and  physiology  of  bone  is  given  by  McLean  (1958).  Herrman 
and  Liebowitz  (1972)  and  Kummer  (1972)  have  written  comprehensive  review 
articles  on  the  mechanical  properties  of  living  bone. 

First  we  will  discuss  the  gross  anatomy  and  histology  of  bone  (pertinent 
to  bone  remodeling)  and  then  consider  the  physiological  mechanisms  within 
bone  which  bring  about  bone  remodeling.  In  Subsection  4,  the  external  condi- 
tions which  bring  about  remodeling  and  the  structural  changes  by  which  re- 
modeling manifests  itself  will  be  discussed. 

I.  2 THE  STRUCTURE  OF  BONE 
a.  Gross  Anatomy 

There  are  two  major  histological  classes  of  bone  tissue  which  signi- 
ficantly contribute  to  the  structural  strength  of  the  skeletal  system.  They  are 
called  cancellous  and  cortical  bone.  First,  cancellous  bone  will  be  described. 
Figure  1 shows  a partial  longitudii  al  section  of  a thigh  bone,  or  femur.  To- 
* ward  the  central  axis,  in  the  epiphyseal  region  (near  the  joints),  there  exists 

a network  of  hard,  interconnected  filaments  called  "trabeculae.  " Interspersed 
among  the  trabeculae  are  marrow  (blood  generating  tissue)  and  a large  number 
of  small  blood  vessels.  The  tissue,  consisting  of  the  trabeculae,  the  marrow, 
and  the  blood  vessels,  is  called  collectively  "cancellous  bone,  " "trabecular 
bone, " or  "spongy  bone.  " 

Lining  the  inner  surface  of  the  marrow  cavities  is  a thin  layer  of 
specialized  tissue  called  endosteum.  Similarly,  the  periosteum  covers  the 
entire  outer  surface  of  the  femur,  except  for  a portion  cf  the  femoral  neck. 

The  endosteum  and  periosteum  are  involved  in  the  regeneration  of  bone  when 
a fracture  occurs.  Immediately  below  the  periosteum,  on  the  outer  surface 
of  the  epiphyseal  region,  there  is  a layer  of  harder,  denser  tissue  called 
"cortical  bone.  " Cortical  bone  contains  no  marrow  and  its  blood  vessels 
are  microscopically  small.  The  details  of  its  microstructure  will  be  dis- 
cussed in  the  next  subsection.  As  we  can  see  in  Figure  1,  as  one  moves 
from  the  epiphyseal  region  into  the  diaphyseal  region  (away  from  the  joints), 

(1)  the  trabeculae  of  the  cancellous  bone  become  scarcer  and  thinner  and 

(2)  the  shaft  diameter  of  the  femur  decreases,  while  (3)  the  layer  of  cortical 
bone  of  the  outer  surface  of  the  femur  becomes  thicker.  Thus  we  can  con- 
clude that  the  cancellous  bone  determines  the  structural  properties  of  the 
femur  in  the  epiphyseal  regions  while  the  cortical  bone  determines  the 

• structural  properties  of  the  femur  in  the  diaphyseal  region. 
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Figure  1.  A Longitudinal  Section  of  a Human  Femur  (redrawn 
after  Kraus  (1969)). 
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I.  3 THE  PHYSIOLOGY  OF  BONE  REMODELING 


Remodeling  manifests  itself  as  a change  in  the  structural  configuration 
of  living  bone.  By  necessity  this  change  of  configuration  involves  a chemical 
reaction  brought  about  on  the  microscopic  level  by  the  activity  of  certain  cells 
embedded  within  the  matrix  structure. 

a.  The  Bone  Cells 

Bone  remodeling  can  only  take  place  by  the  activity  of  certain  cells 
embedded  within  its  matrix  structure.  There  are  three  types  of  bone  cells 
which  actively  participate  in  bone  remodeling:  osteoblasts,  osteocytes,  and 
osteoclasts.  We  will  begin  by  discussing  osteoblasts.  Osteoblasts  synthesize 
bone  material.  In  the  first  step  01  this  synthesis,  the  osteoblast  secretes  a 
soluble  protein  called  tropocollagen  into  the  interstitial  fluid.  The  tropo- 
collagen  is  then  polymerized  into  fibers  of  bone  collagen.  Next,  the  collagen 
matrix  is  mineralized  by  the  precipitation  of  hydroxyapatite  crystals  into  the 
collagen  fibers.  According  to  Frost  (1969)  this  process  takes  about  two  weeks. 

When  an  osteoblast  completes  its  building  process,  it  envelopes  it- 
self within  a hollowing,  or  lacuna,  leaving  only  a thin  cytoplasmic  process  in 
contact  with  the  blood  supply.  Internal  changes  then  occur  within  the  osteo- 
blast and  it  becomes  an  ostecyte.  The  main  function  of  osteocytes  is  to 
facilitate  the  transfer  of  calcium  and  other  ions  between  the  bone  mineral 
and  the  blood  plasma  (McLean,  1958).  The  osteocytes  allow  the  haversian 
system  of  bone  to  act  as  a body-wide  ion  exchange  column  helping  to  stabilize 
the  concentration  of  calcium  in  the  blood  plasma. 

Neuman  and  Neuman  (1958)  have  suggested  that  as  an  osteon  ages, 
its  hydroxyapatite  crystals  become  dehydrated  and  lose  their  ability  to  inter- 
act with  the  plasma.  Thus,  to  maintain  the  physiological  function  of  bone 
an  old  dehydrated  osteon  must  be  resorbed,  even  though  its  loading  remains 
constant  and  its  anatomical  structure  remains  intact.  An  osteon  is  normally 
resorbed  by  large  multinucleated  cells  called  osteoclasts.  It  is  thought  that 
osteoclasts  arise  in  vivo  from  the  merger  of  a number  of  osteocytes  or  osteo- 
blasts. Such  a merger  can  be  observed  in  vitro  by  exposing  osteocytes  or 
osteoblasts  to  parathyroid  hormone.  Chemically,  bone  resorption  is  not 
bone  synthesis  in  reverse;  the  collagen  and  the  bone  mineral  are  dissolved 
simultaneously,  the  protein  is  hydrolyzed  and  the  mineral  is  acted  on  by  a 
chelating  agent  and  possibly  a locally  high  hydrogen  ion  concentration. 

b.  The  Control  of  Remodeling 

This  discussion  shows  that  there  are  a number  of  chemical  re- 
actions involved  in  bone  remodeling.  We  turn  now  to  the  question  of  what 
mechanism  controls  these  chemical  reactions.  There  is  no  definitive  answer 
to  this  question;  however,  there  are  two  phenomena  that  are  likely  candidates: 
one  is  electrical  and  the  other  is  chemical. 


b.  The  Microstructure  of  Bone 

Histologically,  all  bone  tissue  can  be  classified  as  connective  tis- 
sue (Maximow  and  Bloom,  1957).  One  of  the  distinctive  features  of  connective 
tissue  is  thepresence  of  an  abundant  extracellular  material  or  matrix.  In  bone, 
not  only  does  the  matrix  account  for  virtually  all  the  mechanical  strength  of 
the  bone  tissue,  but  also  the  volume  fraction  of  the  matrix  is  considerably 
larger  than  the  volume  fraction  of  the  bone  cells.  In  the  next  section  we  will 
discuss  the  activity  and  control  of  the  bone  cells,  and  the  manner  in  which 
they  affect  remodeling  by  modifying  the  microstructure  of  the  bone  matrix  in 
which  they  are  embedded. 

The  primary  histological  unit  of  bone  structure  is  the  osteon.  In 
cortical  bone,  the  osteons  have  a cylindrical  shape  and  average  150  micra  in 
diameter  and  one  to  two  centimeters  in  length.  These  osteons  have  axes, 
most  of  which  run  parallel  to  the  axis  of  the  bone  itself,  and  these  osteons 
are  cemented  together  into  bundles  by  interstitial  lamellae.  Each  osteon  has 
fluid-filled  hollowing  or  lumen  along  the  center  of  its  long  axis.  The  diameter 
of  the  lumen  is  approximately  40  micra.  Within  each  lumen  there  is  a blood 
vessel  that  supplies  nourishment  to  the  bone  cells  in  the  vicinity  of  that  osteon. 
These  blood  vessels  and  lumina  are  all  interconnected  into  what  is  called  a 
haver sian  system.  During  the  development  of  an  osteon,  bone  generating 
cells,  or  osteoblasts,  located  in  the  canal  deposit  a sequence  of  layers  of 
lamellar  bone  on  the  inner  wall  of  the  osteon,  making  the  canal  successively 
narrower. 

Lamellar  bone  is  a composite  material  similar  in  structure  to  fiber- 
glass. There  is  a strong  but  brittle  inorganic  fraction  (crystals  of  hydroxy- 
apatite) reinforced  by  a matrix  of  a tough  but  flexible  organic  fraction  (mostly 
the  protein  collagen).  The  hydroxyapatite,  like  the  glass  filaments  in  the 
fiberglass,  is  oriented  in  a locally  uniform  direction  or  grain.  The  collagen, 
unlike  the  amorphous  rtsin  of  fiberglass,  is  also  composed  of  filaments. 

The  collagen  fibers  run  parallel  to  the  long  axis  of  the  hydroxyapatite  crystals. 

In  cancellous  bone  the  osteons  compose  the  trabeculae.  These 
osteons  do  not  have  any  uniform  shape.  They  are  often  fenestrated  plates  and 
sometimes  cylinders  or  truncated  cones,  connected  together  into  a three- 
dimensional  network  called  the  trabecular  structure.  A typical  trabecular 
structure  is  illustrated  in  Figure  Z.  The  osteons  of  cancellous  bone  are,  in 
general,  not  hollow.  The  bone  cells  are  nourished  by  blood  vessels  outside 
the  trabeculae  running  within  the  trabecular  structure.  Bone  resorption  and 
reconstruction  can  thus  be  expected  to  take  place  from  the  outside  of  the  os- 
teons rather  than  from  the  inside. 
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In  1953  Fukada  (1957)  discovered  that  bone  demonstrates  a piezo- 


electric effect.  He  later  went  on  to  demonstrate  that  this  was  due  to  the 
collagen  and  he  subsequently  measured  the  piezoelectric  matrix  of  both  bone 
and  collagen.  Thus,  there  is  at  least  one  way  in  which  strain  in  bone  can 
generate  an  electric  field. 

In  1962  Basset  and  Becker  proposed  a physiological  mechanism 
for  the  remodeling  of  bone  in  terms  of  its  piezoelectric  properties.  There 
are  at  least  two  ways  in  which  the  presence  of  an  electric  field  can  affect 
bone  growth.  Becker  and  Murray  (1970)  have  recently  reported  that  an  elec- 
tric field  is  capable  of  activating  the  protein  synthesizing  organelles  in  osteo- 
genic cells  of  frogs.  The  presence  of  an  electric  field  near  polymerizing 
tropocollagen  will  cause  the  fibers  to  orient  themselves  perpendicular  to  the 
lines  of  force.  Basset  and  Pawlick  (1964)  have  reported  that  if  a metal  plate 
is  implanted  adjacent  to  a living  bone,  a negative  charge  on  the  plate  will 
induce  the  deposition  of  new  bone  material  on  this  negative  charged  electrode. 
However,  they  reported  that  there  was  no  resorption  near  the  positive  elec- 
• trode.  Thus,  an  electric  field  may  stimulate  the  activity  of  osteblasts  but 

have  no  effect  on  osteoclasts  other  than  to  inhibit  their  activity. 

The  second  candidate  for  a controlling  mechanism  is  chemical,  the 
calcium  concentration.  Calcium  has  been  implicated  as  a controlling  factor 
in  many  other  biochemical  reactions,  for  instance  muscular  contraction. 

Justus  and  Luft  (1970)  have  demonstrated  that  straining  bone  will  increase 
the  calcium  concentration  in  the  interstitial  fluid.  He  showed  that  this  was 
• due  to  a change  in  the  solubility  of  the  hydroxyapatite  crystals  and  also  dem- 

Ionstrated  this  mechanochemical  effect  with  inorganic  crystals  in  vitro.  A 

chemical  mechanism  would  be  particularly  convenient  in  explaining  the  stimu- 
lation of  osteoclastic  activity.  We  mentioned  above  that  there  are  two  differ- 
ent situations  which  would  call  for  bone  resorption  and  osteoclastic  activity: 
low  ambient  loading  and  dehydration  of  bone  mineral.  Both  of  these  situations, 
however,  would  have  the  same  chemical  effect:  a lower  calcium  ion  con- 
centration. Thus,  if  osteoclastic  activity  were  activated  by  a low  calcium 
concentration  both  of  these  effects  could  be  corrected  by  this  same  chemical 
mechanism. 

I.  4 THE  PROCESS  OF  BONE  REMODELING 
a.  General  Observations 


We  will  follow  the  distinction  made  by  Front  (1964)  between  sur- 
face and  internal  remodeling.  Surface  remodeling  refers  to  the  resorption 
or  deposition  of  bone  material  on  the  surfaces  of  the  endosteum  and  the  peri- 
osteum. Surface  remodeling  can  be  observed  correcting  an  improperly  set 
fracture  in  the  bones  of  young  children.  For  instance,  if  the  two  ends  of  a 
healing  femur  are  set  slightly  skew  to  each  other,  the  bone,  when  it  is  healed, 
will  tend  to  bow  when  a load  is  placed  on  it.  However,  surface  remodeling 


grow  on  the  concave  surface,  and  excess  bone  material  to  be  eaten  away  on 
the  convex  surface,  so  as  to  correct  the  bowing  effect. 

Internal  remodeling  refers  to  the  resorption  or  deposition  of  bone 
material  between  the  surfaces  of  the  periosteum  and  the  endosteum.  Because 
of  the  rigidity  of  adjacent  matrix  material  between  the  periosteal  and  endosteal 
surfaces,  internal  remodeling  (in  contrast  to  surface  remodeling)  cannot  change 
the  physical  dimensions  of  a bone.  However,  there  is  experimental  evidence 
indicating  that  internal  remodeling  can  change  the  porosity,  ash  content,  x-ray 
opacity  and  total  weight  of  a living  bone,  as  well  as  modify  its  mecahnical 
properties  such  as  modulus  of  elasticity  and  compressive  strength. 

In  this  paper,  we  will  limit  discussion  to  the  internal  remodeling 
of  cancellous  bone.  A mathematical  model  of  the  external  remodeling  of 
bone  has  been  suggested  by  Gjelsvik  (1973).  The  problem  of  internal  remodel- 
ing of  cortical  bone  has  already  been  discussed  by  Hegedus  (1976b),  Hegedus 
and  Cowin  (1 973). 

b.  Remodeling  in  Cancellous  Bone 

The  general  concept  of  stress-or  strain-induced  remodeling  was 
first  introduced  by  the  German  anatomist  Julius  Wolff  (1892)  and  is  often 
called  Wolff's  Law.  Wolff  pointed  out  that  the  thickness  and  orientation  of 
the  trabecular  plates  of  cancellous  bone  seemed  to  optimize  its  weight  bearing 
capacity  for  the  particular  load  configuration  to  which  it  was  normally  sub- 
jected. 

One"  of  the  best  examples  of  this  load-optimized  structural  pattern 
can  be  seen  in  a radiograph  of  the  proximal  femur,  whose  loading  con- 
figuration is  indicated  in  Figure  3.  Koch  (1917)  analyzed  these  trabecular 
patterns  and  concluded  that  they  could  be  categorized  into  tensile  and  com- 
pressive groups  which,  under  ambient  loading  conditions,  carry  tensile  and 
compressive  loads,  respectively.  These  trabecular  groups  and  their  ana- 
tomical names  are  indicated  in  Figure  4. 

The  largest  load  component  applied  to  the  human  femur  is  the 
downward  compression  applied  by  the  hip  to  the  femoral  head.  In  the  femoral 
head,  this  compressive  load  is  carried  primarily  by  the  Principal  Compres- 
sive Group  of  trabeculae  that  extends  from  the  proximal  end  of  medial  dia- 
physial shaft  to  the  superior  surface  of  the  head.  Since  the  loading  from  the 
hip  is  not  colinear  with  the  axis  of  the  femur  itself,  there  exists  a bending 
moment  in  the  femoral  neck  that  tends  to  turn  it  medially.  This  bending  is 
opposed  by  the  Principal  and  Secondary  Tensile  Groups  of  trabeculae  that 
extend  from  the  proximal  end  of  the  lateral  diaphysis,  through  the  neck,  to 
the  medial  head.  The  Secondary  Tensile  Group  of  trabeculae  receives  auxil- 
iary support  from  the  Secondary  Compressive  Group  that  arises  from  the 
medial  diaphysial  shaft.  There  exists  a fifth  trabecular  group  of  lesser 
structural  importance  in  the  vicinity  of  the  Greater  Trochanter.  This 
Greater  Trochanter  Group  carries  the  vertical  and  medial  forces  from  the 
muscles  which  are  attached  to  the  Greater  Trochanter.  In  Section  V we  will 
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Compressive 


Figure  4.  Diagramatic  Representation  of  the  Five  Normal  Groups 
of  Trabeculae  of  the  Upper  Femur  (redrawn  from  Singh 
et  al.  , 1970). 
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demonstrate  a model  which  will  develop  these  patterns  when  it  is  subjected  to 
the  loading  configuration  of  a human  femur. 


c.  Immobilization  Experiments 

One  of  the  earliest  quantitative  descriptions  of  internal  remodeling 
in  humans  was  conducted  by  Dietrick  et  al  (1948).  Human  volunteers  were 
immobilized  from  the  waist  down  in  plaster  casts  for  periods  of  from  6 to  8 
weeks.  During  this  study  their  urine,  feces,  and  blood  were  analyzed  for 
certain  organic  components  such  as  creatinine  as  well  as  certain  inorganic 
components  such  as  calcium  and  phosphorus.  Four  days  after  the  plaster 
casts  were  removed,  they  resumed  normal  activity.  The  chemical  analysis 
indicated  that  during  the  immobilization  their  bodies  suffered  a net  loss  of  the 
bone  materials  calcium  and  phosphorus.  After  normal  activity  had  been  re- 
sumed, the  mineral  loss  phenomenon  was  reversed  and  the  body  regained 
calcium  and  phosphorus.  It  took  approximately  si-x  weeks  to  recover  from  the 
effects  of  this  immobilization.  When  the  recovery  period  had  ended,  the 
excretion  levels  of  calcium  and  phosphorus,  had  decreased  to  values  indica- 
ting that  the  body  had  reattained  homeostasis  with  respect  to  these  minerals. 

A similar  net  loss  of  calcium  was  reported  by  Mack  et  al  (1967)  for 
astronauts  subjected  to  weightlessness  for  periods  of  four  days  or  more.  Mack 
correlated  this  calcium  loss  with  a decrease  in  the  x-ray  opacity  of  the  bones. 

It  has  been  pointed  out  by  Brannon  et  al.(1963)  that  the  skeletal  loss  in  recum- 
bency or  weightlessness  could  be  partially  reversed  if  the  subject  performed 
certain  exercises  that  apply  loading  to  the  bones,  which  compensates  for  the 
* loading  that  would  have  occurred  under  conditions  of  normal  gravity.  Mack 

has  suggested  that  the  exercises,  included  in  the  regimen  of  the  Gemini  VII 
flight,  were  responsible  for  the  fact  that  the  astronauts  of  this  14-day  mission 
actually  suffered  less  bone  loss  than  the  astronauts  of  the  8-day  Gemini  V 
mission. 

Immobilization  studies  have  also  been  conducted  on  experimental 
animals.  Kazarian  and  von  Gierke  (1967)  reported  an  immobilization  study 
on  rhesus  monkeys  similar  to  that  of  Dietrick.  The  monkeys  were  divided 
into  two  groups.  One  group  was  immobilized  in  whole-body  casts,  while  a 
control  group,  given  an  identical  diet,  was  allowed  to  have  normal  activity. 

After  60  days  both  groups  were  sacrificed.  An  autopsy  indicated  that  the 
immobilized  animals  were  suffering  from  osteoporosis,  a pathological  con- 
dition of  severely  weakened  bone.  When  bone  samples  from  each  group  were 
excised,  it  was  observed  that  the  bone  from  the  nonimmobilized  animals 
was  stiffer  and  stronger,  and  had  a higher  radiographic  density  than  the  bone 
from  the  constrained  animals.  They  concluded  that  the  immobilized  animals 
suffered  a massive  resorption  of  the  trabeculae  of  their  cancellous  bone  due 
to  the  lack  of  strain  stimuli. 
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d.  The  Effects  of  Cyclic  Loading 

Another  significant  qualitative  study  of  bone  remodeling  was  inspired 
by  the  theoretical  work  of  Pauwels  (1948,  1950).  Pauwels  started  with  the 
premise  that  the  diaphysis  of  the  human  femur,  in  normal  bipedal  activity 
such  as  walking  or  standing,  will  be  subjected  to  a nonuniform  load  distribution. 
As  mentioned  earlier,  the  central  axis  of  the  femur  is  not  directly  over  the 
center  of  mass  of  the  body  and  must  be  pressed  medially  (inwardly)  in  order 
to  prevent  the  leg  from  buckling  outward.  Pauwels  theorized  that  this  medial 
bending  moment  is  partially  offset  on  the  femur  by  a constant  lateral  (outward) 
bending  moment  caused  by  the  iliotibial  tract.  Consequently,  a human  femur 
is  normally  subjected  to  a periodically  medial  and  lateral  bending  moment, 
superimposed  on  a periodical  axial  compression.  This  net  loading  configura- 
tion is  illustrated  in  Figure  5.  Pauwels  concluded  that  the  adaptive  properties 
of  cortical  bone  would  cause  the  medial  and  lateral  side  of  the  femur  to  be- 
come stiffer  and  stronger  than  the  anterior  and  posterior  sides  (the  front  and 
back).  Thusthe  cortical  bone  cf  the  human  femur  modifies  the  distribution 
of  its  material  strength  in  a manner  that  corresponds  to  its  load  distribution. 
This  distribution  is  somewhat  analogous  tothat  of  an  I-beam. 

Amtmann  (1971),  in  a study  on  autopsied  human  femora,  verified 
such  a distribution  of  material  properties.  The  cortical  bone  on  the  medial 
and  lateral  sides  of  the  femur  was  observed  to  be  stiffer,  stronger,  less 
porous  and  more  radio-opaque  than  cortical  bone  from  the  posterior  and 
anterior  portions.  This  cyclic  loading  effect  points  out  an  important  property 
of  bone  remodeling,  namely,  that  bone  accretion  is  not  a linear  function  of  the 
load  applied  but  rather  it  is  some  nonlinear  function  such  that  cyclic  loading 
causes  faster  bone  accretion  than  static  loading  of  the  same  average  magni- 
tude. 

e.  Senile  Osteoporosis 

As  a person  ages,  the  ability  of  the  body  to  synthesize  collagen 
markedly  decreases.  Since  collagen  is  an  essential  component  in  the  syn- 
thesis of  new  bone,  aging  slows  down  bone  accretion.  If  bone  resorption  con- 
tinues unabated,  bone  will  become  weaker  and  more  compliant,  even  though 
the  loading  remains  fixed.  This  disease  is  called  senile  osteoporois  and  is 
an  important  cause  of  debilitation,  particularly  in  elderly  women 

I.  5 CONCLUSIONS 

Living  bone,  in  addition  to  its  structural  utility,  also  fulfills  a vital 
physiological  function.  It  acts  as  a storage  depot  for  minerals,  in  particular 
calcium  and  phosphorus.  Since  the  body  is  constantly  using  calcium  and 
phosphorus,  bone  is  constantly  being  resorbed.  Under  normal  loading  condi- 
tions, the  excessive  removal  of  bone  material  in  a particular  location  will 
induce  a greater  than  normal  strain  around  that  point,  which  will  initiate  the 
accretion  of  new  bone  material.  This  constant  remodeling  of  bone  has  a 
beneficial  effect  for  the  structural  utility.  It  optimizes  the  distribution  of 
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FRONT  VIEW 


EFFECTIVE  LOAD 
ON  DIAPHYSIS 


Figure  5.  The  Effect  of  Cyclic  Loading  on  the  Diaphysis  of  a Right 

Femur  during  Walking  (partially  redrawn  after  E.  Gardner 
et  al. , 1969 ). 
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bone  material  both  generally  and  locally.  Generally,  if  bones  are  too  large 
or  too  high  in  calcium  content  (and  thus  too  dense),  there  will  be  excessive 
weight  and  the  organism  will  suffer  a competitive  disadvantage  from  the  loss 
of  mobility.  If,  on  the  other  hand,  the  bone  is  allowed  to  be  too  small  or 
fragile,  there  is  a danger  of  fracture,  in  which  case  the  organism  would 
suffer  an  even  more  precarious  loss  of  mobility.  Locally,  this  optimization 
assures  that  the  various  regions  of  a bone  undergoing  roughly  the  same 
average  strain  carry  their  proportionate  share  of  the  load. 

This  remodeling  process  is  not  fool-proof.  An  astronaut  in  space,  with- 
out normal  gravity,  may  suffer  a massive  loss  of  bone  material  due  to  a lack 
of  normal  stress  on  his  bones.  This  condition  becomes  dangerous  if  he  returns 
to  earth  and  suddenly  attempts  to  resume  normal  activity.  Even  elderly  peo- 
ple on  earth  can  suffer  massive  bone  resorption  due  to  a malfunction  of  the 
normal  bone  remodeling  mechanism.  Development  of  a mathematical  model 
of  the  stress  adaptation  of  bone  would  be  helpful  to  both  groups.  For  the 
elderly  patient,  a mathematical  model  would  allow  a better  understanding  of 
the  etiology  of  senile  osteoporosis.  For  the  astronaut,  it  would  give  medical 
personnel  a diagnostic  tool  to  predict  at  what  point  bone  atrophy  becomes 
i dangerous.  It  might  also  be  used  as  a therapeutic  tool  to  determine  the  types 

of  exercise  that  would  be  most  useful  in  correcting  bone  loss. 


18 


SECTION  II 

THE  BASIC  GOVERNING  EQUATIONS  AND  THE  MODEL 
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II.  1 INTRODUCTION 

A mathematical  model  for  the  stress-adaption  of  cancellous  bone  has 
been  developed  in  which  cancellous  bone  is  represented  as  a chemically 
reacting  porous  medium.  The  solid  constituent  of  the  porous  medium  repre- 
sents the  matrix  structure  of  the  cancellous  bone  and  the  perfusant  constituent 
represents  the  bone  cells,  the  blood  vessels,  and  the  interstitial  fluid.  The 
strain-adapting  properties  of  bone,  caused  by  the  bone  cells,  is  represented 
by  strain-controlled  chemical  reactions  that  transfer  mass,  momentum,  energy, 
and  entropy  between  the  solid  constituent  and  the  perfusant. 

We  will  begin  our  analysis  by  developing  a continuum  theory  of  an 
N-component  solid  mixture  similar  to  those  formulated  by  Truesdell  (1960)  and 
later  by  Bowen  (1968).  The  next  step  will  be  to  set  the  concentration  of  each  of 
these  N "chemical  components"  equivalent  to  the  thickness  of  an  infinitesimal 
trabecular  plate  in  an  arbitrary  direction  and  consider  the  limit  as  N goes  to 
infinity.  The  resulting  General  Linear  Governing  Equations  are  similar  to 
those  characterizing  classical  linear  elasticity  except  that;  (1)  there  exists  a new 
taisor  field  representing  the  trabecular  development  whose  rate  of  change  is  a 
function  of  the  strain;  and  (2)  the  compliance  in  the  stress  - strain  relationship 
is  a function  of  the  trabecular  development. 

In  the  next  section  these  generalized  governing  equations  will  be  special- 
ized to  the  case  of  two-dimensional  adaptation  to  piane  stress.  The  conditions 
of  two-dimensional  adaptation  will  be  relatively  easy  to  simulate  in-vivo,  and 
will  allow  the  collection  of  experimental  data  which  will  be  essential  in  quanti- 
tive analysis  of  bone  remodeling. 


II.  2 THE  DESCRIPTION  OF  THE  MODEL 

This  section  presents  the  basic  rationale  for  the  mathematical  model  of 
bone  remodeling.  It  begins  by  summarizing  certain  basic  biological  and 
physical  observations  of  the  process  of  bone  remodeling.  It  then  describes  how 
the  salient  features  of  these  processes  can  be  incorporated  into  a mathematical 
model  for  the  stress-adaptation  of  living  bone. 

We  summarize  the  salient  physical  and  biological  properties  as  follows: 

(1. ) There  are  three  basic  components  of  whole  bone:  the  bone  cells,  the 

extracellular  fluid  and  the  solid  extracellular  material,  which  is  called 
the  bone  matrix. 

(2.)  The  bone  matrix  is  a solid  structure  with  interconnected  pores.  The 
mechanical  properties  of  whole  bone  are  essentially  the  same  as  the 
mechanical  properties  of  the  matrix. 
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(3.)  The  trabecular  development  of  the  bone  matrix  is  affected  by  the  ambient 
long-term  strain  history  of  the  bane. 

(4.)  The  trabecular  development  of  the  bone  matrix  is  changed  by  the 

addition  of  mass  or  the  removal  of  mass  from  the  bone  matrix.  This 
transfer  of  mass  occurs  as  the  result  of  a chemical  reaction  which  is 
mediated  by  the  bone  cells.  Momentum,  energy  and  entropy  may  also 
be  transferred  to  or  from  the  matrix  structure  by  these  chemical 
reactions. 

(5.)  The  rates  of  these  chemical  reactions  depend  upon  the  strain  and  are 
very  slow.  It  can  be  estimated  from  the  studies  of  Frost  (1969),  using 
tetracycline  labeling,  that  the  characteristic  time  for  such  a reaction 
normally  is  in  the  order  of  months. 

(6.)  The  extracellular  fluid  is  always  in  contact  with  the  blood  plasma.  The 

plasma  supplies  the  materials  necessary  for  the  synthesis  of  bone  matrix. 

Considering  these  observations,  we  will  assume  the  load -adapting 
properties  of  living  bone  can  be  modeled  by  a chemically  reacting  porous 
medium  in  which  the  rate  of  reaction  is  strain-controlled.  The  porous  medium 
has  two  constituents:  ( 1 ) a porous  elastic  solid  representing  the  matrix 
structure  of  bone,  and  (2)  a perfusant,  which  represents  the  cells,  extracellular 
fluid  and  the  blood  plasma  that  flow  through  the  matrix  structure.  A schematic 
diagram  of  this  model  is  shown  in  Figure  6.  The  fact  that  living  bone  is  encased 
in  a living  organism  is  reflected  in  the  model  by  setting  the  porous  structure  in 
a bath  of  the  perfusant.  When  necessary,  we  will  assume  the  perfusant  bath  to 
be  an  isothermal  heat  reservoir,  an  assumption  that  appears  to  be  easily 
justified  by  common  knowledge  concerning  living  organisms.  The  mechanical 
load  is  applied  directly  to  the  porous  structure  across  the  walls  of  the  perfusant 
bath  as  illustrated  in  Figure  6.  The  system  consisting  of  the  porous  structure 
and  its  perfusant  bath  is  considered  to  be  closed  with  respect  to  mass,  heat 
energy,  and  entropy  transfer,  but  open  with  respect  to  momentum  transfer  from 
loading  and  its  associated  energy  transfer.  The  system  consisting  of  only  the 
porous  structure  without  its  entrained  perfusant  is  open  with  respect  to  momen- 
tum transfer  as  well  as  mass,  energy,  and  entropy  transfer.  We  take  the  bone 
matrix  as  our  control  system,  since  the  mechanical  properties  of  the  bone 
matrix  alone  determine  the  mechanical  properties  of  bone.  We  shall  therefore 
write  balance  and  constitutive  equations  for  only  the  bone  matrix.  The  perfusant 
will  be  accounted  for  only  insofar  as  it  transfers  mass,  momentum,  energy,  or 
entropy  to  the  bone  matrix.  The  rate  at  which  these  transfers  of  mass,  momen- 
tum, energy,  and  entropy  occur  will  depend  on  the  local  strain  and  other  indepen- 
dent variables. 

The  model  contains  a further  simplification.  As  the  porosity  of  the  porous 
structure  changes,  the  area  of  the  interface  between  the  porous  structure  and  the 
perfusant  will  also,  in  general,  change.  There  is  not  a direct  relation  between 
the  porosity  and  the  area  of  the  interface.  We  shall  consider  here  only  porosity 
changes  and  we  shall  not  introduce  the  area  of  the  interface  into  the  model  as  a 
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* 


eat  and  mass  seals 


I 


variable.  Also,  we  shall  not  consider  the  thermodynamic  properties  of 
interior  points  in  the  matrix  structure.  The  thermodynamics  of  the  interface 
in  a chemically  reacting  porous  material  has  been  discussed  in  a recent  paper 
by  Bazant  (1972)  on  concrete  curing. 

We  shall  also  make  two  assumptions  concerning  the  momentum  balance 
of  the  porous  structure.  First,  we  shall  assume  that  the  force  of  interaction 
between  the  porous  structure  and  the  perfusant  is  small  compared  to  the  load 
applied  to  the  porous  structure.  The  second  assumption  will  be  that  only 
mechanically  quasi-static  processes  are  considered.  This  is  justified  because 
the  characteristic  time  of  remodeling  is  so  large  compared  to  a characteristic 
time  for  inertia  effects. 

II.  3 THE  NOTATION 

In  this  subsection  we  will  introduce  the  notation  employed  in  our  develop- 
ment of  continuum  mixture  theory.  In  subsection  4 we  will  develop  balance 
equations  for  mass,  momentum,  and  energy  and  an  entropy  inequality,  and  in 
subsection  5,  constitutive  constraints  for  an  N component  solid  mixture  will 
be  obtained.  In  this  development  we  will  consider  balance  equations  and 
governing  equations  that  apply  only  to  the  matrix  without  the  perfusant.  The 
effect  of  the  internal  perfusant  is  accounted  for  by  transfer  terms  in  each  of 
the  balance  equations. 

We  will  begin  by  postulating  the  existence  of  a field  p indicating  the  bulk 
mass  density  of  the  matrix.  The  field  p has  values  such  that  for  any  arbitrary 
region,  R,  the  mass  within  that  region,  M.^,  will  always  be  given  by  the 
integral 


M = f P dv 

R ft 


(1) 


where  (iv  is  the  element  of  volume.  The  solid  matrix  is  assumed  to  consist 
of  N distinct  chemical  components.  With  each  of  these  components  we  will 
associate  a component  bulk  mass  density,  pa,  such  that  within  an  arbitrary 
region  R,  the  mass  of  the  a-th  component  will  always  be  given  by  the  integral. 


mr  = $ P*  dv 

R R 
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and  thus 


T « „ (4) 

i-  p = p • 

a =1 

Next  we  will  consider  the  kinematic  notation.  The  position  of  the 

particle  X® , of  the  a -th  component  at  time,  t,  will  be  described  by  the 
K 

functional  relation 

X?=Xf(XK'  t)  (5) 

where  i = 1,  2,  3 and  K = I,  II,  III  are  indices  labeling  the  three  different 
coordinate  directions  in  the  deformed  and  reference  configurations,  respec- 
tively. Throughout  this  entire  development  we  will  assume  that  the  N com- 
ponents of  a matrix  do  not  move  with  respect  to  one  another.  Thus,  in  all 
subsequent  development  we  will  consider  a single  function  x^  - X^(X,  t)  as 
describing  the  motion  of  all  components  of  the  mixture.  The  velocity  of  the 
particles  is  given  by 


and  the  deformation  gradients  by 


iK  bXr 


In  index  notation,  a field  quantity  followed  by  a subscripted  comma  and 
an  index,  indicates  a gradient  with  respect  to  the  material  coordinates.  Thus 
the  velocity  gradient 


L = grad  v 


can  be  written  as 


L..  = v.  . 
ij  i.  3 


(8) 


A superimposed  dot  indicates  the  material  or  substantial  time  derivative. 
For  any  arbitrary  tensor  quantity  f we  will  have 


9f 

<v< 

— — + V • grad  f 

at  ~ & ~ 


or,  in  index  notation, 


f. 

l 


af. 

r~  + v.  f.  . 
dt  J i,J 


(9) 


Note  that  in  (9)  we  have  employed  the  summation  convention,  that  is,  that 
over  product  terms  with  repeated  subscript  indices  summation  is  to  be  per- 
formed. This  summation  convention  however  will  rot  apply  to  the  superscript 
indices  indicating  the  chemical  component.  In  this  report  we  will  denote  this 
latter  index  by  Greek  letters. 

Applying  (9)  to  the  definition  (7),  we  can  see  that  the  time  derivative 
of  the  velocity  gradient  will  have  the  form 


= L..F.__ 

i J JK 


(10) 


It  can  also  be  shown  that 


div  v det  F = det  F 


(11) 
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or,  in  index  notation. 


v,  , det  F = L,  , det  F = det  F , 
k,  k ~ kk  ~ ~ 


(12) 


where  det  F is  the  determinant  of  the  matrix  of  deformation  gradients  F 

iK 

The  symmetric  part  of  the  velocity  gradients  is  called  the  rate  of  deformation 

tensor  and  is  denoted  by  D .. 

U 


D.. 

U 


= T (L..  + 
2 U 


V 


(13) 


The  above  notation,  as  well  as  that  which  follows,  is  generally  that 
employed  by  Truesdell  and  Toupin  (I960),  Truesdell  and  Noll  (1965)  and 
Bowen  (1970),  with  the  important  exception  that  we  employ  Greek  superscripts 
to  indicate  the  index  of  the  chemical  component  while  the  previous  authors  have 
used  Gothic  subscripts. 

II.  4 THE  BALANCE  EQUATIONS  AND  THE  ENTROPY  INEQUALITY 
a.  Mass  Balance 

Let  us  now  postulate  that  the  conservation  of  mass  equation  for 
the  jj-th  constituent  takes  the  form 


I Padv  = -/vpVds+  J2«  dv 
9t  R 8R  R 


(14) 


where  3R  is  the  surface  of  R,  ds  is  an  element  of  R»  and  n.  is  a unit  normal 
of  9R  and  ctt  is  the  rate  at  which  constituent  a is  being  generated  by  the 
chemical  reaction.  If  we  apply  the  divergence  theorem  and  (9),  (14)  can 
be  written  in  the  form 


a . 0. 

O ■f  p V 

K i,  i 


= c 


a 


(15] 
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The  conservation  of  mass  for  the  matrix  as  a whole  can  be  written  to  be 


I J 


L 

| 


« 


B_ 

9t 


/ p dv 
R 


BR 


v.p  n.  ds  + 
1 1 


/ c dv 
R 


(16) 


where  c1  is  the  total  rate  at  which  matter  is  being  generated  from  all  chemi- 
cal reactions.  Again  the  divergence  theorem  and  (9)  give  the  field  equation 


p + p v.  . = c (17) 

i,  * 

If  we  substitute  (15)  into  (17)  and  (4),  we  will  have 

N 

c = r (is) 

<2=1 


Finally  we  note  that  it  will  be  useful  for  our  future  development  if  we 
define  two  new  sets  of  field  variables  which  can  be  used  in  lieu  of  p and  p °. 
These  new  fields  will  be  called  the  zero  strain  bulk  matrix  mass  density  and 
the  zero  strain  bulk  component  mass  density,  defined  as 


and 


y = p det  F 


a a 

V = P 


det  F , 


(19) 


respectively.  In  terms  of  these  new  field  variables,  the  conservation  of  mass 
equations  (17)  and  (15)  can  be  written  as 


y = c det  F 


.a 

Y 


J f IT 

c det  F , 


respectively. 
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(20) 


b.  Momentum  Balance 


The  conservation  of  momentum  equation  for  the  a-th  constituent  is 
given  by  the  relation 


— f p v.  dv  = -£  (p  v.)  v.n.  ds  + £ T. . n.  ds  + f pttb^  dv 

9t,/ri  -r  i ii  ^iiJ  J v i 


9B 


B 


A®  a“ 

+ f (P  .+  c v.)  dv  (21) 

^ 1 1 
B 


where  T denotes  the  partial  stress  tensor  of  thea-th  component,  b.  is  the 

body  force  per  unit  mass  of  the  a-th  component,  and  p^  is  the  force  applied 

on  the  a_th  component  by  all  other  components  and  the  perfusant.  The  first 

three  integrals  on  the  right  hand  side  of  (8)  are  customary,  the  fourth  is  a 

departure.  The  first  term  in  the  fourth  integrand,  p?;  represents  the  force 

applied  on  the  a_th  component  by  the  (N-l)  other  components  and  the  per- 

aCK 

fusant;  the  second  term,  c v^,  represents  the  momentum  added  to  the  porous 
matrix  structure  due  to  the  addition  of  mass  by  the  chemical  reaction. 

Using  the  divergence  theorem,  the  conservation  of  mass  (15), 


where  T is  the  total  stress  tensor  for  the  matrix  constituent,  defined  as 


a 

T..  = E T..  (24) 

lJ  a=l 


b.  is  the  total  body  force  on  the  matrix,  defined  as 

\ - Q (25) 

and  p.  is  the  force  applied  to  the  matrix  by  the  perfusant.  Again  use  of  the 
divergence  theorem  and  the  conservation  of  mass  equation  (4)  yields  the  field 
equation 


p x.  = T. . , + p b.  + p 
1 1J.  J i i 


(26) 


Comparison  of  (22)  with  (26)  gives  us 


N 

A „ aCK 

P£  = S P,  (27) 

a = i 


Again,  in  the  special  case  when  A*  = 0 for  all  a,  equation  (26)  reduces 
to  the  familiar  conservation  of  momentum  equation  for  a simple  continuum. 

c.  Moment  of  Momentum  Balance 

We  will  next  consider  balance  of  moment  of  momentum.  In  terms  of 
direct  notation,  balance  of  moment  of  momentum  for  the  O-th  component  can 
be  postulated  to  be 


dT~J  ^ * (p 

R 


v )]  dv  = - / (x  x p v)  (v  ds)  + / X * [t  • ds] 

8R  ~ 8R~  ~ ~ 

+ * <P  £ + cv  + p)+m]dv 


(28) 
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where  ds  is  an  element  of  the  surface  of  BR  times  a unit  normal;  the  binary 
operand  x denotes  the  cross  product  of  two  vectors;  and  m is  the  momentum 
supply  vector  for  the  a-th  component  playing  a role  in  (28)  analogous  to  the 

role  of  pi  , in  equation  (2  1).  Using  the  divergence  theorem,  and  equations  (15) 
and  (22),  equation  (28)  reduces  to  the  form 


r a,t 
[T  ] 


= M 

'V 


or 


a 

a 

a 

T. . 

- T.. 

= M.. 

U 

J1 

u 

(29) 


where  the  components  of  the  skew- symetric  tensor 


a 

M.. 

U 


are  given  by 


and 


M 


M 


M 


[?, 

= m“  = 

1 1 

22 

a 

a 

= -M 

23 

32 

a 

= - M** 

31 

13 

a 

d 

= - M , 

12 

21 

33 
= m. 


0! 

*3 


(30) 


We  will  now  postulate  a conservation  of  moment  of  momentum  equation 
for  the  matrix  as  a whole  to  be 


J x (pv)  dv  = - ^ ( p v x v)(v  • ds)  + fa  x x [T  • ds] 
R~  ~ 9R~  ~~6R~ 


+ _/*  x (p  b + p*cv)  ) + m]  dv  (31) 

where  m is  the  momentum  supply  vector  for  the  matrix.  Equation  (18)  can 
be  reduced  to  the  field  equation 


A 


(32) 


T..  - T..  = M.. 
1J  Ji  iJ 


where 


a N 

M..=  l 

1J  cT=l 


M 


a 

ij 


(33) 


Of  course,  if  all  the  moment  of  momentum  supply  vectors  are  zero,  equation 
(32)  reduces  to 


T..  = T .. 
Ji 


(34) 


which  is  the  conservation  of  moment  of  momentum  equation  for  a simple 
material. 

d.  Energy  Balance 

In  this  subsection  we  consider  conservation  of  energy.  The  conserva- 
tion of  energy  equation  of  the  a-th  component  has  the  form 


(35; 


OC  Ot 

where  £ is  the  specific  internal  energy  of  the  a-th  component;  q^  is  the 

heat  flux  vector  in  the  a-th  component,  rtt  is  the  specific  heat  supply  for  the 

— 0( 

a-th  component  per  unit  time  and  h is  the  energy  transfer  between  the  matrix 
and  the  (N-l)  other  components  and  the  perfusant.  The  first  three  integrals 
on  the  right  side  of  (35)  are  customary:  the  first  integral  is  the  sum  of  the 
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convected  internal  and  kinetic  energy;  the  second  is  the  work  done  by  the 
surface  tractions  and  the  heat  flux;  and  the  third  represents  the  heat  supply 
per  unit  mass  plus  the  work  of  the  body  force.  The  fourth  integrand  in  (35) 
is  associated  with  the  energy  transfer  to  the  a-th  component  from  the  (N-l) 
other  components  and  the  perfusant.  The  first  term  in  the  fourth  integrand  is 
the  work  done  by  the  force  of  interaction;  the  second  represents  the  kinetic 
energy  of  the  added  mass;  and  the  third  represents  the  internal  energy  of  the 
added  mass,  assuming  its  specific  internal  energy  is  the  same  as  that  in  the 


extant  porous  structure.  The  fourth  term  represents  any  energy  transfer 

between  the  a-th  component  which  is  not  accounted  for  by  the  other  three 

— Ci 

terms.  For  example,  h may  represent  the  energy  transfer  due  to  the  fact 
that  the  internal  energy  of  the  added  mass  is  not  the  same  as  that  of  the  ex- 
tant material. 


Using  the  divergence  theorem,  and  equations  (15)  and  (26),  we 
can  rewrite  (35)  as 


a ;a  = Ta  La  - qa  . + ^ ra  + ha 

u i. 1 


We  will  next  postulate  the  energy  balance  equation  for  the  matrix 


as  a whole  to.be 


TT  / P (e  + J v v ) dv  = - / p (e  + j v v )v  n ds  + / (v  T - q ) n ds 

c i 1 OD  ^ 1 i 1 J 1 ‘J  J J 


+ J p (bv.  + r)  dv  + J (p.v.  + ~£~  + ce  + h ) dv 


where  we  have  defined 


t 


N 

= (i/p)  z p 


a a 

e 


a=l 


N 

qi = 2 
a=i 


a 

qi 


N 

and  r = (l/p)  Z Pttra  (38) 

a=i 


We  should  observe  that  (38)  is  made  possible  by  the  fact  that  the  N chemical 
components  of  the  matrix  are  not  moving  relative  to  one  another.  We  can  next 
use  the  divergence  theorem  and  the  equations  (17)  and  (26)  to  yield  the  field 
equation 


p e = T..L,..  - q.  . + p r + h 
ij  J1  i. 1 


where 


- y -a  y a A 
h=lh  + Zce  - c e 


O'  =1 


a =i 


(39) 


(40) 


e.  The  Entropy  Inequality 


With  each  component  O.  we  will  associate  a component- specific  entropy 
function  rj  , such  that  within  any  region  R the  entropy  of  the  a-th  chemical 
component  is  given  by 


(41) 
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Similarly,  for  the  matrix  as  a whole,  we  will  assume  that  there  exists  a 
matrix- specific  density  function  rj  such  that  with  any  region  R the  entropy 
of  the  matrix  is  given  by 


SR  = / Pt!  dv  = t S®  • 

r a =1 


(42) 


These  two  assumptions  place  no  constraints  on  the  generality  of  this 

development.  Next  we  will  postulate  an  entropy  inequality  for  each  component 

such  that  the  partial  time  derivative  of  each  S has  the  form 

R 

8 r a a , . 

iT  J p n dv  - 


f 

8R 


paqa  v.n.  ds 


1 l 


/ 

8R 


r a _ 

q.  n. 

1 1 

ds  + / 

P r 

a.,  x f 

r Ta  1 

* at  a , h_ 

L « 

a 

dv  + J 

C q + a 

0 

R 

L e J 

R 

L e J 

dv 


(43) 


a 

where  0 is  a positive  real  scalar,  indicating  the  temperature  of  the  a-th 

= 01  — Ql 

component  and  h is  defined  in  a manner  analogous  to  that  of  h in  equation 

= a —a 

(35).  We  distinquish  between  h and  h to  indicate  that  all  the  energy 
transferred  to  the  CK-th  chemical  component  need  not  contribute  to  the  entropy 
production  in  the  at -th  component.  As  in  the  case  of  the  balance  equations  for 
energy  and  entropy,  the  first  three  integrals  on  the  right  side  are  traditional 
and  the  fourth  integral  is  a departure:  The  first  integral  represents  the  con- 
vected  entropy,  the  second  is  the  entropy  production  associated  with  the  heat 
flux,  the  third  is  the  entropy  production  associated  with  the  specific  heat 
supply  and  the  fourth  represents  the  entropy  production  due  to  heat  transfer 
and  chemical  reactions  with  the  (N-l)  other  chemical  components  and  the 
perfusant.  The  first  term  in  the  fourth  integrand  is  the  entropy  produced 
by  the  chemical  reaction  generating  new  material  in  the  a-th  component  as- 
suming that  this  new  material  has  the  same  specific  entropy  as  the  extant 
material.  The  second  term  in  the  fourth  integrand  represents  all  entropy 
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production  not  already  accounted  for.  For  example,  this  term  may  represent 
the  additional  entropy  contribution  due  to  the  fact  that  the  temperature  of  the 
newly  formed  matrix  material  has  a higher  temperature  than  that  of  the 
extant  material. 


Equation  (43)  may  next  be  converted  into  a field  equation  by  the  use  of 
the  divergence  theorem,  conservation  of  mass  (15)  and  (22)  giving  us 


V 

1 

- ■« 

a ot  , ra 

o r + h 

a 

+ 

a 

0 

0 

, i 

(44) 


Using  the  same  motivation  as  that  above,  we  can  postulate  that  the 
entropy  inequality  for  the  matrix  as  a whole  has  the  form 


JL  f 

at  j pT1 

R 


N 


dv  > ^ pq  v.n.  ds  - £ Yj 

0R  11  9R  a = l 


q®  n." 

l l 


L 0 


Ot  J 


N 

*•*11 

r a=i 


a a A a .a-  -a 

p r + cri  9 + h 


(45) 


which  reduces  to  the  field  equation 


N 


Ph  > 


a=l 


u*\ 

OL  Ot  , 7 OL 

p r + n 

a 

+ a 

\ e I 

0 

, i 

J 

(46) 


In  further  discussions,  we  will  limit  ourselves  to  the  special  case  when 
all  N chemical  components  have  the  same  temperature  field  0 , in  which  case 
inequality  (45)  can  be  written  in  the  form 


Ph  > 


where 


h=  £ =h*+  l ?Ve  -c> 

a-i  a = i 


(47) 


(48) 


dv 
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However,  this  inequality  does  not  contain  any  new  restrictions  that  are  not 
already  contained  in  inequalities  (44). 

f.  The  Helmholz  Free  Energy 

In  this  final  subsection  we  will  recast  the  conservation  of  energy  and 
the  entropy  inequality  in  terms  of  the  Helmholz  free  energy.  Let  us  define 
the  c omponent- specific  free  energy  density  as 


ol  a.  at  a 

4j  = £ - T]  0 


and  the  matrix -specific  free  energy  as 


n Ol  C 

4>  = ( l /p  ) X p 4* 

a =1 


Since  we  have  assumed  that  the  temperature  fields  of  all  N chemical  com- 
ponents of  the  matrix  are  the  same,  we  can  also  write 


= e — r|  © 


Finally  we  can  substitute  equation  (47)  into  the  energy  balance  equation  (36) 
and  obtain  the  relation 


a • at  a • . a „ . a ot  -a  a a 

-p  M + q 9 + r)  6 ) - T L + q - h = p r 

J1  1 


for  conservation  of  energy  in  the  OS-th  component.  In  addition,  equation  (52) 
can  be  substituted  into  (44)  to  give  the  entropy  inequality 


a • a a -ot  „ jot  T * at  /l\  at  ^ n 

"P  + - p ^ 6 + Tij  Lji  + h - (ej  qi  gi  - 0 


L.L ... 


II.  5 THE  CONSTITUTIVE  CONSTRAINTS  AND  THE  ENTROPY  INEQUALITY 


a.  Introduction 

In  the  previous  subsection  we  discussed  the  balance  equations  for  mass, 
momentum,  energy,  and  entropy.  Now  we  make  constitutive  assumptions 
appropriate  for  an  elastic  material  and  obtain  constraints  on  the  constitutive 
equations  from  the  entropy  inequality  and  from  certain  physical  and  physiolog- 
ical assumptions  about  the  physiological  nature  of  bone  remodeling.  The 
resulting  equations  constitute  the  governing  equations  for  an  N-component, 
ideal,  "non-swelling",  linearly  elastic'  solid  mixture.  In  the  following  sub- 
section we  will  discuss  how  this  finite  element  component  model  might  be 
generalized  to  simulate  cancellous  bone. 

In  the  following  subsection  we  will  summarize  the  balance  equations 
and  the  constitutive  constraint  that  can  be  built  into  the  balance  equations.  In 
the  next  subsection  we  will  discuss  the  manner  in  which  the  entropy  inequality 
placed  further  constraints  on  the  constitutive  equations.  In  the  fourth  sub- 
section we  will  discuss  certain  auxiliary  constitutive  assumptions,  based  on 
simple  physical  and  physiological  assumptions  about  the  nature  of  bone  and 
bone  remodeling,  which  will  further  simplify  the  governing  equations. 
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b.  The  Balance  Equations 

We  will  recall  that  in  our  development  of  the  balance  equations,  it 
simplified  our  discussion  to  place  certain  constraints  on  the  governing  equa- 
tions. These  constraints  were,  that 

(1)  all  N components  of  the  matrix  have  the  same  motion,  i.  e. 


a 


x (X,  t)=  x(X,  t)  for  a = 1,  2,  ...  N 


(55) 


(2)  all  N components  of  the  matrix  have  the  same  temperature 

6=  6P( X,  t ) = 0 (X,  t)  for  a=  1,  2,  . . . N (56) 

It  would  also  be  physically  reasonable  to  expect  further  that 

(3)  all  N components  of  the  matrix  have  the  same  body  force,  i.e. 


ra 


b = b (X,  t ) = b(X,t)  for  tt=  1,  2 . . . N 
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(57) 


.!  i 


By  letting 


N 

Y = Zy* 


T..=  £ T® 
1J  a=l  1J 


N 

* -z  < 

' a = ! 


v a ( 

Yr  = h Y r 

a =1 


N 

h=i>“ 

a=i 


V ^ = Z Y0^** 

a =1 


y*i  = Z V*ha  i 

-1  J 

the  conservation  of  mass,  momentum,  and  energy,  for  the  matrix  as  a whole 
takes  the  form 


Y = c det  F 


v(v  - b ) = T det  F 
i i ij  ~ 


T - T . = 0 

ij  j1 


and 


(I  + q0  + T,  e - r)  ± (T  L ~ \i+  h>  det  Z > 0 


As  mentioned  previously,  given  the  constraint  (56),  the  entropy  inequality 
for  the  matrix  as  a whole  does  not  contain  any  new  constraints  not  already 
implied  by  (60). 

c.  The  Entropy  Inequality 

Our  objective  is  to  model  the  normal  adaptive  processes  that  occur  in 
cancellous  bone  remodeling  as  strain-controlled  mass  deposition  and/or 
resorption  which  modify  the  trabecular  structure  of  the  porous  elastic  solid. 
We  shall  therefore  make  constitutive  assumptions  similar  to  those  made  for 
elastic  solids,  but  with  the  inclusion  of  additional  independent  variables 
|y^  , yZ  , . . . , y^  } which  measure  the  degree  of  trabecular  adaptation.  Thus 
we  can  write 

zy  AQf  O 

* =4-  (6,  e .,  f iK , yp) 

Ti°  = ^a(e  . e(i’  f.k,  yP) 


Ta.  = Ta.  (6  , e . , F , /) 
1J  1J  » 1 iK 


= e t , F.K,  /) 


a OL  (-iOf  /a  rt  -r- » 

c = C (9  * 6 , i ’ FiK’  7 > 


ha  = ha(e  , e . . f.k,  y&) 


a at  AQt  R 

h = h*(0  , 0(.  , F.k,  yp) 


X 2 n 

Where  the  third  argument  of  (63)  symbolizes  the  set  y , y , y 

Substituting  (63)  into  (6o^  gives  us  the  entropy  inequality 


“ Ul  e + ^ 


86 


a . Dt  Nr 

e + y 

86  . ,i 

,i  a-1 


2^  c*  det  f"1  l..f  +„<V 

8/  ~ 8FiK  )K  ! 


[a  Aa  /i\  ft 

T . . L ..  + h - {-)  q . 6 . det  F > 0 

ij  ji  \6/  -1  .1  J — 


(64) 


where  equations  (62),  and  (10)  have  also  been  employed. 


We  will  now  use  an  argument,  initially  outlined  by  Coleman  and 
Noll  (1963)  and  Coleman  and  Mlzel  (1964);  we  will  claim  that  the  time  de- 
rivative terms  6 


and  L..  can  be  specified  independently  of  the  in- 
XJ 

o 0 o 

dependent  variables  (0,  0 F , yp ).  Note  however  that  yP  is  constrained 
by  equations  (60)^  and  (63)’^.  The  entropy  inequality  (64)  then  implies  the 
following  restrictions 


24*  + a = 0 , = 0 , and  t"  = 

86  + 11  ’ 86  . 13  det  F 


9FiK 


FjK 


(65) 


These  results  show  that  ^“is  independent  of  the  temperature  gradient  and  the 
t,®  and  Ttt  are  related  to  iJj®  by  the  formulae 

ij 


Aft 

--  4*  ( e 


F y ^ ) 
*iK’7  ' 


9 


a 


and 


8^tt 


8F . 


iK 


(66) 
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When  the  results  (66)  are  substituted  back  into  (65),  the  reduced  entropy  in- 
equality becomes 


Sir  + h“ 
8/ 


(67) 


From  equation  (61)^  it  follows  that  the  total  stress  for  the  matrix  material  has 
the  form 


T.. 

U 


a 


34^ 


3F. 


iK 


) 


(68) 


This  result  will  prove  to  be  important  in  future  discussion. 

d.  Approximations  and  Additional  Constraints  on  the  Governing  Equations 

The  governing  equations  derived  in  the  previous  subsubsection  are  too 
general  to  have  immediate  practical  application  and  it  will  be  necessary  to 
• simplify  and  specialize  them  further.  These  simplifications  are  justified  by 

certain  easily  observable  physiological  properties  of  living  bone.  For 
instance  we  can  assume  that 

(1)  living  bone  may  be  regarded  as  isothermal,  i.  e. 

0 = e“(X,  t)  = constant  for  or  = 1 , 2, . . . N 

and  all  X and  t . (69) 


Under  these  circumstances,  the  governing  equations  of  the  previous  sub- 
section reduce  to  the  set 


y*  = cU  (F.k,  yP)  det  F 


and 


y(v.  - b.)  = T..  . det  F , 
' i i iJ,J  ~ 


U 
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where 


N 

T a 

r*  L > . 

a = l 


We  will  next  assume  that 

( 2 ) bone  remodeling  occurs  in  a manner  such  that  bone  remains 
nonswelling. 

By  nonswelling  we  mean  that  as  trabecular  adaptation 
proceeds,  the  geometric  configuration  of  the  bone  under  zero  stress  (the  re- 
ference state  for  strain)  remains  the  same,  no  matter  where  the  trabecular 
adaptation  takes  place.  One  might  imagine  a block  of  cancellous  bone  con- 
taining four  equidistant  points,  the  vertices  of  a tetrahedron  marked  for  the 
purpose  of  measuring  strain.  When  adaptation  occurs,  material  is  added  or 
taken  away  from  the  trabeculae,  but  if  the  boundary  conditions  are  ever 
returned  to  conditions  such  that  the  local  stress  within  the  strain-tetrahedron 
is  zero,  the  distance  between  the  four  vertices  of  the  tetrahedron  will  always 
return  to  their  initial  value.  In  terms  of  the  stress  function  in  eq.  (63)  this 
constraint  can  be  written  as 

0 = Ta.  (6  y ^ ) for  all  7/3  (71) 

ij  iK 


and  in  terms  of  the  free  energy  function  we  will  have 


°-W~.  *“(FiK-  y“)  '• 


iK 


F =6 
iK  iK 


(72) 
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Let  us  next  define  the  strain  tensor  by 


EKM  2 EFiKFiM  ‘ 6KM^  (73) 

In  terms  of  the  displacement  gradients  u^  ^ ERM  has  the  representation 

EKM  “ 2 tUK,  M ' UM,  K " UP,  MUP,  (?4) 

In  addition,  the  free  energy  function  (63)  and  the  rate  of  reaction  function 
in  eq.  (63)  can  be  redefined  in  terms  of  these  new  arguments  as 

_a  B />0£  B 

c (Ekm-  > > - c (FiK'  yp' 

(EKM'  yP)  = (FiK’  ] (75) 


respectively. 


We  are  now  in  a position  to  make  the  assumption  that 

(3)  bone  is  subjected  only  to  small  strain  and  this  strain  is 
linearly  related  to  the  stress. 

The  constitutive  equations  for  c and  0 can  now  be 
approximated  by  assuming  that  the  displacement  gradients  ^ are  small 

and  their  squares  can  be  neglected.  When  the  squares  of  the  displacement 
gradients  are  neglected  and  the  material  and  spatial  coordinate  systems  are 
made  to  coincide,  the  strain  tensor  can  be  approximated  as 


E..  = 
ij 


I [ u.  . + u.  . ] , 

L LJ  J,ij 


(76) 


a ~a 

Approximations  for  c and  0 are  obtained  by  neglecting  terms  of  order 
three  and  higher  in  a Taylor  series  expansion 


_a 

c 


(E  ) - aa(y^)  + A“  (y^)  E..  + 1/2  B*  (y^)  E. . E, 

U U ' ij  ijk/  ' 1 tj  k/ 


U 
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*°'V  yf)  * *“ (0’ yf>)  + Dl  + Kjki  <*“  > Ea Ek« 


where 


a.  /3  . _a  . £ 

a (r  ) = c (y  , 0) 


a , j3,  d a fi  , 

Ay  ^ '“rET^  (EPq-  ? > i 

IJ 

E =0 

pq 


a B d a R 

Bn*y > = ra-  *■  <v  ^ > 

ij  k£ 


E =0 

pq 


Dijiyl>)  = Je-  *“(e  ./)  ! 

ij  J 


E =0 

pq 


E =0 

pq 


The  constraint  (72)  requires  that 


a B B 

D. . (yH  ) = 0 for  all  a and  y 

ij  T 

thus  when  (79)  is  substituted  back  into  (77),  the  stress-strain 
relationship  (70)  takes  the  form 


1 


/s 

T..  = C.. 


ij  Cijkjt 


The  approximations  of  the  current  subsection  allow  us  to 
simplify  considerably  the  governing  equations  developed  in  the  previous 
subsection.  The  approximated  governing  equations  consist  of 

(1)  a definition  of  the  strain  tensor  in  terms  of  the  displacement, 


E..  = | [u.  . + u.  .]  ; 
ij  L i.J  h 


(2)  a stress-strain  relationship 

Ty  - S*.  > Ek« 


where 


N 

a R ~ a a 

cijki  <>  > = 2,  * cijk^ ; 

oc=l 


(82) 


(3)  a rate  of  reaction  equation  for  each  component 


»0f 


a 


ij 


va 

1J  1J 


r (E;;,  yH)  = a“  + a“  e.  . + ±b“  , e.  . e. 


ijk/  ij  k l 


where 


where 


(d)  a conservation  of  momentum  equation 


via.  = T. . . + yb.  , 
a aj.J  ' a 


y-L 


In  order  to  quantify  the  physical  properties  of  this  material  it  is 

necessary  to  specify  two  sets  of  constitutive  functions,  each  of  which  is 

specified  by  its  own  set  of  material  parameters.  The  first  set,  composed 

solely  of  the  function  , describes  the  total  stiffness  of  the  material 

as  a function  of  the  chemical  composition  of  the  material.  For  this  reason 

the  parameters  C?.,  will  be  called  the  Stiffness  Parameters  and  the  function 

fi..,  . will  be  calleuthe  Stiffness  Function.  The  second  set  of  functions  ca 
ijk£ 

describes  the  rate  at  which  the  chemical  composition  of  the  material  changes 
as  a function  of  the  strain  and  the  current  chemical  composition.  Since 
these  changes  in  the  chemical  composition  are  brought  about  by  the  remodel- 
ing  properties  of  bone,  aQ  , A--  , and  might  be  termed  Remodeling 

Parameters  and  the  function  set  ca  may  be  said  to  collectively  describe  a 
Remodeling  Function. 

These  equations  are  similar  to  those  of  a classical,  linear  elastic 
material  except  that  there  now  exists  an  N parameter  set  {ya}  describing 
the  chemical  composition  of  the  elastic  material  whose  values  are  deter- 
mined by  the  conservation  of  mass  equations  (82).  In  fact,  in  the  special 
case  when  the  Remodeling  Function  is  identically  zero,  i.  e.  , ca  = 0 for  all 
tt  , then  the  chemical  composition,  as  parameterized  by  the  parameter  set 
{ ya}  remains  fixed  and  the  governing  equations  reduce  to  those  of  a 
classical  linear  elastic  material.  In  the  case  of  living  bone,  however,  the 
Remodeling  Function  is  not  zero  and  this  chemical  change  in  the  composition 
allows  the  material  to  modify  or  adapt  its  physical  properties  to  its  ambient 
loading  conditions.  Consequently,  the  governing  equations  (82)  may  be  said 
to  describe  an  Adaptive  Elastic  Material. 


II.  6 THE  CANCELLOUS  BONE  MODEL 
a.  Introduction 

In  the  previous  subsection  we  developed  an  adaptive  elastic  model  in 
which  there  were  N "chemical  components".  In  this  subsection,  we  will 
introduce  a model  for  cancellous  bone  in  which  these  "chemical  components" 
are  rendered  equivalent  to  the  trabecular  plates  that  can  be  observed  and 
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measured  anatomically.  We  will  start  off  in  the  second  subsubsecticn  by 
equating  these  N "chemical  components"  to  N infinitesimal  trabecular 
plates  oriented  normal  to  N trabecular  basis  vectors.  Since  we  do  not  know 
a priori  the  spatial  orientation  of  the  trabecular  basis  vectors,  we  will 
reformulate  the  governing  equations  so  that  the  dependence  of  the  N spatially 
oriented  chemical  components  is  replaced  by  a trabecular  density  function 
whose  values  will  be  independent  of  the  N trabecular  basis  vectors. 

In  the  third  subsubsection  we  will  demonstrate  that  this  trabecular 
density  function  can  be  approximately  specified  by  a symmetric,  second - 
rank  tensor,  M , called  the  trabecular  adaptation  tensor.  The  tensor  field 
M has  a very  convenient  physical  interpretation:  Anatomically,  the  matrix 
material  of  cancellous  bone  is  arranged  into  a structure  of  mutually  ortho- 
gonal trabecular  plates.  The  eigen  vectors  of  M may  be  interpreted  as 
indicating  the  normal  of  each  of  these  trabecular  plates,  while  the  eigen  value 
of  each  of  these  eigen  vectors,  indicates  the  relative  thickness  of  these 
trabecular  plates.  In  Section  IV,  we  will  discuss  design  of  instrumentation 
which  will  enable  investigators  to  measure  the  tensor  M experimentally. 


b.  The  Trabecular  Distribution  Function 

We  will  begin  by  associating  each  of  the  N "chemical  components" 
of  the  previous  subsection  with  a trabecular  plate  whose  normal  is  the  unit 
vector  w0'.  The  set  of  all  N w®  may  be  referred  to  as  the  trabecular 
basis  vectors,  since  they  play  a role  somewhat  analogous  to  the  basis  vectors 
in  the  matrix  representation  of  tensors.  The  components  of  any  trabecular 
basis  vector  in  a three-dimensional  space  can  be  resolved  as 


a a a 

w i = cos  (B  ) sin  (0  ) 

w^  = sin  ( 9 a)  sin  (0**)  (84) 

a a 

W^  = COS  (0  ) 


where  6 and  0 are  the  longitude  and  colatitude  parameters  of  conventional 
spherical  polar  coordinates.  These  angles  in  three-dimensional  space,  with 
relationship  to  the  unit  vector  w®,  are  illustrated  in  Figure  7.  We  will  also 
note,  for  future  reference,  that  the  unit  vector  i_  3 can  be  rotated  into  the  unit 
vector  w®  by  left  multiplying  by  the  orthagonal  transformation 

[cos  0 cos  0j  [-sin  0 cos  0]  [-sin  0] 

£ sin  0 ] [ cos  0 ] [ 0 ] 

[cos  0 sin  0 ] [ sin  0 sin  0]  [ COS  0]| 


Q (0  ,0)  = Q (w)  = 


(85) 
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We  are  now  in  a position  to  introduce  the  trabecular  distribution 
functions  v (w)  and  V (0,  0)  which  can  be  used  in  lieu  of  the  variable  set 
{ ya}.  Let  us  consider  a unit  sphere  covered  with  a thin  shell  of  mass  y 
which  is  composed  of  material  with  a uniform  mass  density  of  y.  We  will 
further  assume  that  the  surface  of  this  unit  sphere  has  been  partitioned  into 
N sectors  each  of  which  is  associated  with  one  of  the  unit  vectors  wa.  Such 
a partitioning  is  illustrated  in  Figure  8.  If  we  consider  the  area  of  each 
sector  to  be  AS®  and  the  mass  of  the  material  in  each  sector  to  be  y®,  then 
the  quotient 


V 


a - 


y 


a. 


As® 


(86) 


can  be  considered  to  be  approximately  proportional  to  the  thickness  of  the 
shell  in  the  sector  a.  We  next  observe  that  the  definition  (86)  will  still 
remain  valid  regardless  of  how  the  sectors  are  partitioned  or  even  of  how  the 
trabecular  basis  vector  set  { w®}  is  selected.  This  suggests  that  we  consider 
the  set  of  all  v ® depending  on  w®  as  defining  a trabecular  distribution 
function  v (w)  such  that 

va  = V ( w)  I (87) 

~ QC 

w = w 

whose  domain  is  the  set  of  all  unit  vectors  yv  . The  trabecular  distribution 
function  can  also  be  expressed  as  a function  of  the  longitude  and  the  colatitude 


VU  = V (0,0)  | 

0 = 0 a 

0 = 0“ 


(88) 


Using  eq.  (86),  the  expression  for  the  total  mass  density  can  be  written  as 


N 

^ _ a a ,nn. 

y = v (w  ) AS  ; (89 ) 

a- 1 


while  the  expression  for  the  total  stiffness,  indicating  the  Stiffness  Function, 
takes  the  form 


/N 

Cijkf 


N 

' a — a . a 
C. ..  . v (w  ) AS  ; 
— ljke  ~ 

a=  1 


(90) 
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and  equations  specifying  the  rate  of  trabecular  adaptation  take  the  form 


_d 
dt 


["(w)]  • = c“(E..;  AS1!/1,  AS2!/2,  AS3!/3, . . . ASNi/N). 


(91) 


As  the  number  of  trabecular  plates,  N , approaches  infinity,  each  segmental 
area  ASa  becomes  incrementally  small  and  the  summations  (89)  and  (90) 
become  the  surface  integral  expressions 


v (w  ) dS 

r^j 


(92) 


and 


CijkX 


rsj 


Hw01)  dS 


(93) 


Thus,  as  N approaches  infinity,  expressions;  (89)  and  (90)  and  (91)  indicate 
that  the  values  of  y and  an<^  ^ (w)  depend  on  every  value  of  the 

function  V (w)  when  the  argument  ^ is  a unit  vector.  In  general,  a function 
whose  value  depends  on  the  val  ues  of  a second  function  over  a given  subdomain 
of  this  second  function  is  called  a functional.  The  simplest  example  of  a 
functional  is  the  conventional  integral.  An  integral  can  be  thought  of  as  a 
function  whose  value  depends  on  a second  function,  the  integrand,  over  a 
subdomain  of  the  integrand  as  indicated  by  the  limits  of  integration.  There- 
fore, given  that  the  total  bulk  mass  density  and  the  total  stiffness  are  functionals 
of  the  trabecular  distribution  function,  it  is  not  surprising  to  see  them  in  (87) 
and  (9  3)  expressible  in  terms  of  weighted  intergrals  of  the  trabecular  distri- 
bution function  for  all  values  of  w such  that  | y/j  = 1. 

However,  since  functionals  may  not  always  be  expressible  as  weighted 
integrals,  we  now  introduce  a new  notation  in  which  a functional  is  indicated 
by  an  oversized  script  letter  under  which  is  placed  an  expression  containing 
a dummy  variable  prescribing  the  domain  of  functional  dependence  of  the 
functional  argument.  For  example,  in  this  notation  the  expression  in  eq.  (9  3) 
can  be  recast  into  the  form 
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This  new  notation  is  necessary  to  describe  the  rate  of  trabecular  adaptation. 
In  the  case  of  expression  (91),  as  N approaches  infinity,  the  time  derivative 
of  v (w)  depends  on  not  only  w and  E but  also  on  all  values  of  the  function 
V (w)  for  | ^/|  = 1.  Thus  in  the  new  notation  the  expression  (91)  approaches 
the  form 

7(w)  = & jw;E;V(w')(  . (9  5) 

~ |V  |=1  1 ~ ' 

as  N approaches  infinity. 


The  next  step  in  our  development  of  the  cancellous  bone  model  is  to 
use  a principle  of  continuum  mechanics  called  material  objectivity.  In 
general,  material  objectivity  requires  that  the  governing  equations  have  the 
same  mathematical  formulation  when  written  with  respect  to  any  arbitrary 
nonaccelerating  coordinate  system.  In  our  case,  material  objectivity  requires 
ihat  equations  (9  3)  and  (94)  will  be  invariant  when  rewritten  in  terms  of  a new 


i • , , ; ..  -i  f ^ 

coorainaie  system  uuidmcu  iium 


the  orthogonal  transformation  (85). 


instance,  if  the  integration  of  expression  (93)  is  performed  with  respect  to  the 
new  coordinate  system,  the  expression  for  the  Stiffness  Function  takes  the  form 


v (w)  Q (w)  Q.  (w)  Q , (w)  Q ,(w)  C 


pqrs 


where  Qjj  (w)  is  given  by  (85)  and  C_  rg  describes  the  stiffness  of  a single 
trabecular  plate  in  the  xi,x2  plane.  Likewise,  material  objectivity  requires 
that  the  functional  expression  (9  5)  take  the  form 


v (w) 


' Q (w)  E Q(w);  v (Q(w ) w7 


|w'|=l 


P 
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where 


{ii’  $ v{™)} 


K 1 = 1 

I r^,  I 


(98) 


and  represents  the  remodeling  function  of  a trabecular  plate  in  the  X2 
plane.  We  might  also  observe  that  the  assumption  of  quadratic  dependence 
of  ya  on  Ejj  called  for  in  eq.  (82)  will  constrain  the  functionalj^p  ° of 
(97)  to  have  the  form 


& |e;  t/(w'  ) J 

I™'  l=1 


|w7  1=1 

i rsj  i 


+ E.. 
ij 


Is' 1=1 


+ E..E 
ij 


ki 


l = 1 


(99) 


where  , dr  f and 


are  functionals  on  the  function 


V (w7 ) over  the  domain  ]w7  J =1  which  are  independent  of  E. 


c.  The  Trabecular  Adaptation  Tensor 

In  the  previous  subsubsection,  we  discussed  how  it  was  possible  to 
specify  the  trabecular  adaptation  of  cancellous  bone  using  the  Trabecular 
Distribution  Function  T7(w).  In  this  subsubsection,  we  will  discuss  how  the 
Trabecular  Distribution  Function  can  be  approximately  specified  by  a 
symmetric,  second  rank  tensor  M called  the  Trabecular  Adaptation  Tensor. 

It  will  be  shown  that  the  Trabecular  Distribution  Function  can  be  approximated 
by  a quadratic  form  V (w)  — sf  M w and  that  the  tensor  M is  subject  to 
physical  interpretation  in  terms  of  the  anatomical  properties  of  the  trabeculae 
of  living  cancellous  bone. 

We  will  begin  by  considering  how  we  might  approximately  specify 
the  Trabecular  Distribution  Function  using  a finite  number  of  parameters. 
Since  the  function  p(0,0)  is  a periodic  function  in  both  its  arguments,  the 
obvious  place  to  look  for  such  approximating  parameters  would  be  in  a double 
Fourier  series.  Let  us  consider  a double  Fourier  series  for  v (6,0) 
truncated  at  terms  of  order  two,  which  has  the  form 
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2 2 


i»(6,  d>)  = ) ) n"111  cos  m6  cos  n 0 + ) ) N™  sin  m8  cos  n0 

' i i cc  L — ' sc 
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V T1  mn  - . . C V mn  . 

> ) N cos  m0  sin  n0  + / ) N sin  m0  sin  rv±> 
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m=0  n=l 


m=l  n=l 


(100) 


where  the  Fourier  coefficients  are  given  by  the  integrals 
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0 ) sin  m0  sin  n0  d0  d0 


(101) 


It  will  also  prove  convenient  to  express  (100)  in  the  form 


V (0,0)  = Aq(0)  + Aj(0)  CO8  0 + A2(0)  cos  20  + B^O)  sin  0 + B2(0)  sin  20 


where 


2 

An(0)  = l (N^n  cos  m®  + N^n  sin  ) 
m=0 


and 


B (0)  = Y (Nmn  cos  m0  + N1^11  sin  m0  ) . (102) 

n L>  \ cs  ss  ' 

m=0 

The  physical  nature  of  the  distribution  of  the  trabeculae  places  mathematical 
constraints  on  the  Trabecular  Distribution  Function  which,  in  turn,  constrains 
the  Fourier  coefficients.  For  instance,  we  can  assume  that  the  Trabecular 
Distribution  Function  is  a continuous  function  of  the  longitude  0 even  at  the 
"north  pole"  and  the  "south  pole"  of  the  sphere  illustrated  in  Figure  7.  This 
requires  that 

^ (0  , n ir)  = constant  for  n = 0,  +1,  + 2,  . . . (103) 


In  addition,  we  can  observe  that  a trabecular  plate,  defined  by  a unit  vector 
w , will  be  indistinguishable  from  a trabecular  plate  defined  by  the  unit  vector 
-w.  This  will  require  that 

P(w)  = V (-w) 

rv  rsj  ' 

and 

A 

v (0  ,0 ) = v (0+tr , 0 + ir)  (104) 


When  these  constraints  are  applied  to  the  series  (102'»  they  take  the  form 
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aq(0)  = Ao(0+ir) 

A (0)  - A (0)  = constant 

u 

Aj  (0)  = 0 

Bjtf)  - 0 

and  B2(0)  = -B2(0+ir).  <105) 

And  when  they  are  applied  to  the  Fourier  coefficients  they  constrain 
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01  10  11  12  22  10  11 
while  they  renuire  the  coefficients  N , N , N , N ,N  ,N  ,N  , 
1 1 cc  cc  cc  cc  cc  sc  sc 

N12,  N21,  N01,  N02,  N11,  N21,  N11,  N12,  N21,  and  N22  to  be  equal  to 
sc  sc  cs  cs  cs  cs  ss  ss  ss  ss 

zero.  This  allows  us  to  rewrite  the  series  (1 00)  into  the  form 


(0 , 0)  = N°°  + N°2  + [N°2  + N22  cos  20  + N22  sin  20  J ^cos  20  - lj 


cc 


sc 


+ |^N  cos  0 + N^2  sin  ©J  J^sin  20  j 


(107) 


where  there  are  only  six  independent  Fourier  coefficients. 


The  next  step  in  our  development  is  to  consider  how  the  function 
V (w)  can  be  expressed  in  terms  of  these  Fourier  coefficients.  If  we 
substitute  the  definition  of  the  vector  components  of  w (84)  into  the 
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truncated  Fourier  series  (107)  , V (w)  can  be  expressed  as  a quadratic 
form 


where 


U (w)  = M. . w.  w. 
~ D i 3 


for  | w I = 1 


(108) 
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M__  = M - N2J 
23  32  ss 


21 

M . = M = N 
13  31  cs 


and 


M _ = = -2N22 

12  21  sc 


(109) 


The  extreme  simplicity  of  the  expression  (108)  suggests  that  we  use  the 
tensor  M to  specify  the  trabecular  adaptation,  in  lieu  of  the  Trabecular 
Distribution  Function.  We  will  call  the  symmetric,  second  rank  tensor  M 
the  Trabecular  Adaptation  Tensor  and  observe  that  it  will  have  the  same 
properties  of  objectivity  and  symmetry  as  the  symmetric  second  rank  tensor 
describing  the  stress  and  strain. 


Furthermore,  like  the  stress  and  strain  tensors,  the  Trabecular 
Adaptation  Tensor  is  subject  to  physical  interpretation.  To  demonstrate 
this  physical  significance,  we  will  first  observe  that,  using  integrals  (100) 
and  constraints  (106),  the  tensor  components  of  M can  be  expressed  directly 
in  terms  of  the  function  17  (w)  by  the  expression 


M.. 

ij 


dS 


(110) 
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Let  us  next  assume  that  the  material  described  by  this  integral  possesses,  at 
any  given  point,  only  a finite  number  of  trabecular  plates.  Using  definitions 
(86)  and  (87)  this  integral  can  then  be  recast  as  the  finite  summation 


Li 


M..  = 
ij 


a a a 

w.  w.  y 
i J 


where  the  integer  N now  refers  to  the  number  of  the  trabecular  plates  at 
the  given  point.  We  next  remark  that  anatomical  observation  of  living  bone 
indicates  that  cancellous  bone  is  generally  arranged  in  a pattern  of  mutually 
perpendicular  trabecular  plates.  Consequently,  if  we  define  an  orthogonal 
tensor  L as  a continuous  function  of  position  such  that  the  coordinate  axes 
are  always  mapped  by  L(x)  into  normals  of  the  trabecular  plates,  we  can 
express  the  Trabecular  Adaptation  Function  in  the  form 

r xi  ° ° l 


M = L 


0 0 L 

2 ~ 


0 OX. 


where  Xi  , X£  an  d X3  describe  the  bulk  mass  density  of  the  trabecular 
plates  normal  to  the  vectors  L _i  j , L and  L i.3  » respectively.  If, 
in  addition,  we  assume  that  the  local  mass  density  of  the  lamellar  bone 
making  up  the  trabeculae  has  a uniform  mass  density,  we  can  make  the  fol- 
lowing physical  interpretation  of  the  tensor  field  M : 

1)  the  eigen  vectors  of  ^1  represent  normals  to  the  trabecular 
plates  of  the  matrix  structure  of  cancellous  bone;  and 

2)  the  eigen  values,  corresponding  to  these  eigen  vectors  are 
proportional  to  the  thickness  of  these  trabecular  plates. 

It  therefore  can  be  concluded  that  the  trabecular  adaptation  of  cancellous 
bone  can  be  described  mathematically  by  evaluating  the  Trabecular  Adaptation 
Tensor  and  this  Trabecular  Adaptation  Tensor  can  be  determined  by  evaluating 
the  fields  Xl  X2  X3  and  L . In  Section  IV  we  will  discuss  procedures  by 
which  measurement  of  these  parameters  might  be  effected. 


d.  The  Constitutive  Functions  and  Material  Objectivity 

At  the  end  of  the  previous  subsection  we  concluded  that  in  order  to 
quantitatively  specify  the  properties  of  an  Adaptive  Elastic  Material  it  would 
be  necessary  to  specify  two  sets  of  functions:  A Stiffness  Function  and  a set 
of  Remodeling  Functions.  In  this  subsection  we  demonstrated  that  it  is 
possible  to  represent  the  trabecular  adaptation,  formerly  specified  by  the 
parameter  set  [ya},  by  a symmetric,  second-rank  Trabecular  Adaptation 
Tensor,  M . In  this  subsubsection  we  will  consider  how  these  Stiffness 
and  Remodeling  Functions  can  be  reformulated  in  terms  of  the  tensor  M 
and  how  the  spectral  decomposition  of  and  material  objectivity  with 
respect  to  M effect  this  reformulation. 

Let  us  first  consider  the  Stiffness  Function,  which  we  wish  to  evaluate 
as  the  function 

C...  . = C...  . (M)  (113) 

ljkf  ljkf  ~ 

If  the  objectivity  conditions  are  applied  to  equation  (96),  the  stiffness  tensor 
takes  the  form 


c..=  y Q.  Q.  Q , Q . F(Qw“)dS 

ilk  A Z_  in  in  rk  a A ~~ 


ijk i L ip  jq  rk  ~si 

ot- 1 


Where  Q is  an  arbitrary  orthogonal  transformation  and  N is  the  (finite) 
number  of  trabecular  plates.  However,  we  have  already  stated  that  we  can 
consider  N=  3 . If  we  now  use  the  approximation  (108)  and  consider  the 
vectors  Qw  as  being  the  eigen  vectors  of  M , then  eq.  (113)  takes  the 
form 


Cijkf 


= C..,  .(M)  = ) X Q.  Q Q , Q . C 

ljkjt  ~ L ip  jq  rk  si  pqrs 
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(116) 


and  where  \a  and  L-  are  given  in  eq.  (119)  and  Cpqrs  describes  the 
stiffness  of  a single  trabecular  plate  in  the  xi  X2  plane. 

Let  us  next  consider  the  equations  describing  the  rate  of  the  chemical 
reactions.  On  rewriting  (95)  in  terms  of  the  variable  M the  resulting 
equation  would  have  the  form 


M 


£ (E,  M) 


(117) 


where  the  second  rank  tensor  function  £ can  be  called  the  Remodeling 
Function.  The  objectivity  requirement,  in  direct  notation,  takes  the  form 


/s 


Q S (QE^  ; QMS)  Q=  &(£;  M) 


(118) 


We  also  know  that  the  strain  dependence  of  6 is  quadratic.  Thus,  analogous 
to  eq.  (99),  eq.  (117)  can  be  recast  into  the  form 


Wij  = 'Mab»  + + EkJEmm  <Mab» 


(119) 


Lf  we  now  apply  the  objectivity  constraint  eq.  (118)  to  this  expression,  we 
then  have 
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M«  = Qim<?  rF  (QcaQdbMcd>  + Qim%°,p°sqEr,£™<°caQdbMcd> 

+ QimQinQrpQsqQv.QwuErsEvw»^u  (QrpQsqM,s» 


for  any  orthogonal  Q.  In  the  special  case  when  Q equals  the  tensor  Q .. 
given  by  the.  formula  in  eq.  (116),  then  the  arguments  of  the  Functions  f'kk 
and  become  diagonalized.  Thus,  (120 ) can  be  recast  again  in 

terms  of  three  new  functions  and  ^mn  which  depends  only  on 

the  eigen  values  of  M,  i.  e. 


M = Q Q7  F^fX  , X , X ) + Q3  Q3  Q3  Q3  E fi””  (X  , X.,  XJ 
lj  1m  jn  12  3 lm  jn  rp  sq  rs  pq  1 2 3 


+ q3  Q3  Q3  Q3  Q3  Q3  E E A™  (X  , X_,  X ) 
Uim  jn  rp  sq  rt  wu  rs  vw  pqtu  1 2 3 


2 1 

Furthermore,  if  we  also  consider  Q and  Q we  will  have  the  constraints 
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e.  Summary  of  the  Cancellous  Bone  Model 

We  began  this  subsection  with  a mathematical  model  for  cancellous 
bone  in  which  the  trabecular  adaptation  was  described  by  a variable  set 
[y &)  indicating  the  bulk  mass  density  of  a finite  number  of  trabecular  plates. 
This  formulation  was  unsatisfactory,  since  we  did  not  know  a priori  the 
orientation  of  the  trabecular  basis  by  which  or  through  which  these  trabecular 
plates  were  ordered.  In  the  second  subsubsection,  we  replaced  dependence 
on  the  variable  set  {y^}  with  dependence  on  a Trabecular  Distribution 
Function  V (w)  which  gave  our  model  an  infinite  trabecular  basis.  Specifi- 
cation of  the  trabecular  adaptation  in  terms  of  the  Trabecular  Distribution 
Function  did  free  us  from  dependence  on  an  arbitrary  fixed  trabecular  basis, 
but  the  resulting  governing  equations  were  very  cumbersome  to  deploy.  In 
the  third  subsection,  we  demonstrated  that  the  Trabecular  Distribution  Function 
could  be  approxima  ted  by  a quadratic  form  v (w)  = M w and  hence  the 
trabecular  adaptation  of  cancellous  bone  can  be  represented  by  a symmetric, 
second-rank  tensor,  M , called  the  Trabecular  Adaptation  Tensor.  This 
tensorial  representation  of  the  trabecular  adaptation  in  terms  of  M is  very 
convenient  mathematically,  since  it  allows  M to  be  subject  to  the  same 
objectivity  and  symmetry  considerations  that  apply  to  the  tensors  fields  E 
and  J.  In  addition,  the  tensor  M has  a very  convenient  physical  interpreta- 
tion: the  eigen  vectors  of  Jvl  represent  normals  to  the  trabecular  plates  while 
the  eigen  values  of  Jyl  represent  the  thickness  of  these  plates. 

The  governing  equations  for  cancellous  bone  consist  of 
a)  a definition  of  the  strain  tensor  in  terms  of  the  displacement, 


E..  = Hu.  - +u.  .]  ; (122) 

ij  2 L i. J J.i 


b)  a stress-strain  relationship 


T.. 

ij 


(123) 


1 


62 


and 


(d)  a conservation  of  momentum  equation 


V u.  = T. . . + y b. 
i XbJ  i 


where 


y = — tr  M ; 
2 ~ 

n 


where  M is  symmetric  and  hence  can  be  resolved  as 


M = iL4 
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with 


L L = 1 


(126) 


(127) 


(128) 


(129) 


To  specify  a traction  boundary  value  problem,  the  following  data  must  be 
given 
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1)  The  geometric  configuration  of  the  body  B and  the  two  con- 
stitutive functions,  the  Stiffness  Function  (124)  and  the  Remodeling  Function 
(125); 


2)  The  time  interval  (0,to)  for  which  a solution  is  desired; 

3)  The  body  force  field,  b,  acting  on  the  body  B; 

4)  The  surface  traction  _s  specified  on  SB  for  (0,tQ),  where 


s - n T (x,  t)  for  x C SB 
i J ij  ~ ~ 


(130) 


and  nj  is  a unit  normal  to  the  surface  SB  ; 

5)  The  continuous  initial  displacement  field  v°(x),  the  continuous 
Initial  Trabecular  Adaptation  Tensor  field  M (x),  and  the  continuous  intital 
velocity  field  v°(jc)  for  all  £ € B where 

~°(~)=!r  u(&  fc)1t-o  (13D 


With  these  data,  the  traction  boundary  value  problem  of  quasi- static,  isothermal, 
adaptive  elasticity  is  to  determine  the  fields  u(x,  t),  M(x,_t),  £(x,  -t),  X(x*-i) 

corresponding  to  b which  satisfy: 

1 ) The  initial  conditions 


£(*/  °)  = 2° <x)>  M(x,  0)  - M°(x)  on  B; 


(132) 


2)  The  boundary  conditions 


n.T..(x.  t)  = s.(x,  t)  on  3B  ; 

j ij^  l ^ 


(133) 


and  equations  (121)  through  (127). 
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II.  7 INTERPRETATION  OF  THE  MODEL 


Let  us  now  interpret  these  equations  in  terms  of  the  observable 
remodeling  properties  of  cancellous  bone.  Assuming  that  bone  in  a particular 
region  initially  has  a trabecular  adaptation  of 


0 

0 


(134) 


when  expressed  in  the  appropriate  coordinate  system.  We  will  first  consider 
the  behavior  of  the  system  under  uniaxial  loading  T^j  = -P°  6^3  6j 3 whose 

magnitude  P°  (>  0 ) is  such  that  the  total  rate  of  remodeling  on  all  trabeculae 
is  zero.  Physiologically,  a zero  rate  of  remodeling  on  all  trabecular  plates 
may  be  termed  homeostasis  and  occurs  when  the  net  rate  a bone  accreation 
on  each  trabecular  plate  exactly  cancels  the  rate  of  bone  resorption.  Mathe- 
matically, homeostasis  requires  that  the  strain  E°  have  a value  such  that 
T^j  = Cijk^(M°)  e£^  and  that  ft(E°;M°)  = 0.  We  will  also  note  in  passing 

that  the  symmetry  will  require  that  X^  = X9,  . 

Let  us  next  consider  what  would  happen  if  the  load  were  increased,  i.  e.  , 
the  value  of  T were  modified  to  T—  = TP.  - £ ^3  . Provided  A>0,  the 

absolute  value  of  the  strain  will  increase  , and,  in  particular,  the  negative 
value  of  the  strain  component  E^  will  become  more  negative.  Physiologically, 
an  increase  in  the  axial  compressive  strain  might  be  expected  to  stimulate 
osteoblastic  activity  in  the  trabecular  plates  which  run  along  the  X3  axis,  i.  e.  , 
the  plates  which  are  perpendicular  to  the  xj  and  x2  axes.  Mathematically,  a 
more  negative  component,  E^,  will  cause  the  Rjj  and  R22  components 
of  the  remodeling  function  to  become  positive  and  the  M^  and  M22  components 
of  the  trabecular  adaptation  tensor  to  have  larger  positive  values.  As  the  two 
trabecular  plates  normal  to  the  Xj,  and  x2  axes  become  thicker,  the  strain, 
induced  by  the  loading  T^j  tends  to  decrease  and  the  net  rate  of  bone  accreation 
tends  to  decrease.  Mathematically,  this  means  that  as  Mjj  and  M22  increase, 
E33  becomes  less  negative  and  the  Rjj  and  R22  components  of  the  remodeling 
function,  which  bring  about  the  trabecular  thickening,  also  decrease.  Thus,  as 
long  as  the  load  T^j  continues  to  be  applied,  the  trabecular  plates  continue  to 
be  thickened,  but  at  a progressively  slower  rate.  Mathematically,  we  can  say 
that  the  system  asymptotically  approaches  a new  state  of  homeostasis  for  the 
ambient  loading  T^j  in  which  the  trabecular  adaptation  tensor  has  the  form 


M = M7 


X7  0 0 

0 X'2  0 

0 0 x3 


(135) 
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where  X*  2 >^<2  » X*  3 > X*3  • In  addition,  symmetry  again  requires 

that  X7j  = X physiologically,  the  attainment  of  this  new  state  of  homeo- 
stasis corresponds  to  an  increase  in  the  axial  loading,  causing  an  increase  in 
the  thickness  of  the  trabeculae  supporting  the  axial  load.  The  effect  of 
increased  loading  is  schmatically  described  in  the  second  frame  of  Figure  9. 

Likewise,  we  can  expect  that  the  system  will  respond  to  a decrease  in 
axial  loading  in  an  analogous  fashion:  Given  a system  where  the  trabecular 
adaptation  tensor  is  initially  M°  given  by  (130)  and  where  the  loading  is 
decreased  to  the  value  Ty  = T?j  - A ^3  ^3  with  A<0,  the  system  can  be 
expected  to  approach  a new  state  of  homeostasis  in  which  the  trabecular  adapta- 
tion tensor  will  have  the  form 


M = M* ' = 
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where  X^  = X/ 3 < X°  = X^  and  X^  < X^  . Physiologically,  this  corres- 
ponds to  a decrease  in  the  axial  loading,  causing  a decrease  in  the  thickness 
of  the  trabeculae  supporting  this  axial  load.  The  effect  of  a decreased  load  on 
the  trabeculae  normal  to  the  X£  axis  is  schematically  illustrated  in  the  third 
frame  of  Figure  9. 


Finally,  we  will  consider  the  response  of  the  system  to  a shear  loading. 
Let  us  consider  the  system,  initially  with  a trabecular  adaptation  M°  given  by 
(130)  , in  which  the  load  is  of  the  form 


T*  = -P 

r^j 


[ sin  0 ] [sin  0 cos  0] 
[ sin  0 cos  0 ] [ cos  0 ] 


(137) 


where  P>0  . It  is  not  difficult  to  demonstrate  that  this  load  is  nothing  more 
than  the  uniaxial  load  T 0 = P 6536:3  rotated  0 radians  around  the  x^  axis. 

If  0 is  small,  cos^0  = 1 and  sin^  0 is  negligible  and  thus  the  principal 
difference  between  the  load  T°  and  the  load  T'1'  is  the  presence  of  the  shear 
component  T^3  and  T32  in  the  latter.  These  shear  stress  components  will, 
in  turn,  cause  the  off-diagonal  components  of  the  remodeling  function  R23 
and  &32  to  be  nonzero.  Physiologically,  the  appearance  of  a shear  load  can 
be  observed  to  induce  remodeling  of  the  trabecular  structure  whose  net  effect 
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Figure  9.  Trabecular  Configuration  Under  Various  Loading  Conditions. 
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is  to  eliminate  the  shear  strain.  This  can  be  accomplished  if  this  remodeling 
caused  reorientation  of  the  trabecular  structure  so  that  the  principal  axis  of 
the  tensor  M is  made  to  coincide  with  the  principal  axis  of  the  ambient  strain 
tensor.  This  means  we  can  expect  that  the  Trabecular  Adaptation  Tensor  will 
asymptotically  approach  the  value 


M*  = M 


X°  0 0 

0 [ X°  cos20  + X°  sin2e]  [(X°  - X°  ) sine  cos6] 


0 [(X°  ~ X°)  sin  0 cos  0 ] [ X°  cos2e  + X°  sin  6] 


(138) 

which  is  nothing  more  than  the  initial  Trabecular  Adaptation  Tensor  rotated 
0 radians  around  the  axis.  The  effect  of  this  shear  loading  is  schemati- 
cally described  in  the  fourth  frame  of  Figure  9. 
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SECTION  III 


THE  PROBLEM  OF  TWO-DIMENSIONAL  STR ESS- ADAPTATION 
AND  ITS  NUMERICAL  APPROXIMATION 

III.  1 INTRODUCTION 

The  previous  section  developed  a mathematical  model  for  cancellous 
bone.  This  sectitn  will  demonstrate  its  application.  However,  any  demonstra- 
tion of  tiiis  model  will  require  that  we  first  must  solve  the  governing  equations 
(122)  - (126  ).  Because  of  the  complexity  of  solving  a system  of  nine  partial 
differential  equations  in  nine  unknowns  (three  displacement  and  six  trabecular 
adaptation  components  ) we  would  be  well  advised  to  seek  plausible  methods 
of  simplifying  their  solution.  In  subsection  III.  2,  we  will  introduce  two 
methods  of  simplification:  the  two-dimensionalization  and  finite  element 
approximation  of  the  governing  equations.  The  former  simplification  is 
motivated  by  the  fact  that  it  is  possible  in  Classical  Elasticity  to  reformulate 
a number  of  technically  important  boundary  value  problems  in  only  two  dimen- 
sions. Finite  element  approximation  allows  us  to  avoid  the  virtually  impossible 
task  of  finding  an  analytic  solution. 

However,  the  use  of  numerical  techniques  has  one  distinct  disadvantage: 

In  order  to  solve  the  simplified  governing  equations  that  are  developed  in 
Subsection  III.  2 we  must  have  explicit  formulations  for  the  analogs  of  the 
constitutive  functions  (124)  and  (125).  Although  there  are  a number  of  sources 
from  which  one  might  plausibly  estimate  the  form  of  the  stiffness  function  (124), 
there  is,  with  a single  exception,  nothing  in  the  literature  from  which  one 
might  determine  the  form  of  the  remodeling  function.  This  single  exception 
is  the  author's  dissertation  (Hegedus  1976),  in  which  a special  case  of  the 
stres s -adaptation  of  cortical  bone  was  analyzed.  In  his  analysis  the  author 
used  certain  well-documented  qualitative  properties  of  bone  remodeling  to  obtain 
phenomenological  constraint  on  the  remodeling  function.  With  these  phenomeno- 
logical constraints  it  was  possible  to  obtain  an  analytic  solution  to  the  problem 
of  the  stres s -adaptation  of  cortical  bone  in  the  special  case  of  uniform  uniaxial 
stress. 

In  subsection  III.  3,  we  will  introduce  finite  element  techniques  and 
use  these  techniques  to  argue  that  cancellous  bone  may  be  approximated  as  a 
pinned  triangular  structure  whose  structural  cross -members  (trabecular  links) 
behave  mechanically  as  miniature  copies  of  the  cortical  bone  analyzed  in  the 
author's  dissertation.  In  subsection  III.  5 we  will  discuss  the  organization 
and  application  of  a series  of  computer  programs  which  simulate  the  stress- 
adaptation  properties  of  cancellous  bone,  using  the  model  subject  to  the  above 
approximations. 
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We  will  show  that  the  model  predicts  that  an  initially  homogeneous 
adaptive  elastic  body,  in  the  shape  of  a human  femur,  subject  to  the  boundary 
traction  analogous  to  that  of  a human  femur  will  develop  values  for  the  tensor 
field  M which  correspond  to  the  trabecular  patterns  which  are  visible  in 
the  human  femur  in  situ. 

HI.  2 PLANE  STRESS 

As  mentioned  above,  in  classical  elasticity  it  is  possible  to  reduce  a 
number  of  technically  important  boundary  value  problems  in  terms  of  governing 
equations  of  only  two  dimensions.  One  of  the  best  examples  of  such  a reduction 
in  dimensions  occurs  in  the  case  of  Plane  Stress  . Plane  Stress  (in  the 
X1  x2  -plane)  is  said  to  occur  if  the  stress  tensor  components  T13  ’ T23  ’ and 

T33  vanish  for  all  x and  t.  The  conditions  of  plane  stress  are  most  easily 
demonstrated  under  the  following  circumstances: 

Let  us  consider  a generalized  cylinder,  as  illustrated  in  Figure  10,  with 
a generator  parallel  to  the  x^  axis  and  with  bases  in  the  planes  x-j  = + h. 

This  body  can  be  said  to  be  'thin'  in  the  x^  direction,  provided  that  the 
height  of  the  cylinder,  2h,  is  small  compared  to  the  linear  dimensions  of  the 
generator.  Let  us  now  assume  that  all  boundary  forces  applied  to  the  surface 
of  the  cylinder  dB  are  applied  only  to  its  lateral  surface  <JB'  and  that  these 
forces  exist  only  in  the  xix2~  plane  and  are  uniform  with  respect  to  the 
variable  x-j  . Let  us  further  assume  that  the  inertial  and  body  lorces 
are  negligible  and  that  the  material  composing  the  cylinder  is  also 
uniform  and  symmetric  with  respect  to  the  x3axis.  Under  these 
circumstances,  if  we  define  the  two-dimensional  "reduced"  stiffness 
tensor 


^abcd  ^abcd 


Cab33  Ccd33 


(139) 


and  the  two-dimensional  "average"  displacement  vector 


x = + h 


u = u (x.  , x ;t ) = — 

a a 1 2 2 


L I u (x,  ,x  ,x  It)  dx, 
2h  / a 1 2 3 : 


(140) 


X3=  -h 


where  the  indices  a,  b,  c and  d take  on  only  the  values  1 and  2 ; then  it  is 
possible  to  show  that  the  governing  equations  (122)  and  (123)  reduce  to  the  set 
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and 


( CabcdEcd>b 


2 E12,12  Ell,  22  + E22, 11 


while  the  boundary  conditions  (133)  can  be  written  in  the  form 


(CabcdEcdhb» 


s (x  ; t) 
a "" 


xedB' 


xedB' 


(141) 


(142) 


Equation  (141)  is  called  the  compatibility  equation  and  assures  the  existence 
of  a vector  field  u"a  such  that 


1_ 

2 


u ) 

b,  a 


(143) 


The  derivation  of  the  compatibility  requirement  is  thoroughly  discussed 
in  Sokolnikoff  (19  56). 

The  governing  equations  (141)_constitute  a set  of  three  equations  in  the 
three  unknowns  En  , E12  > and  E 22  • I11  the  special  case  when  the  boundary 

conditions  (142)  are  independent  of  time,  the  solution  tensor  E will  also  be 
independent  of  time  and  equation  (141)  taken  with  (142)  will  be  said  to  constitute 
a plane  stress  elastostatic  boundary  value  problem.  These  are  a number  of 
techniques  by  which  elastostatic  boundary  value  problems  may  be  solved. 

One  very  elegant  and  powerful  method  using  complex  variables  was  originally 
suggested  by  Kolossoff  and  was  subsequently  developed  by  Muskhelishvili 
(19  53)  . A very  readable  discussion  of  the  use  of  complex  variables  applied 
to  elastostatic  boundary  value  problems  is  found  in  Sokolinikoff  (19  56).  Another, 
even  more  powerful,  albeit  less  elegant,  method  of  solving  the  elastostatic 
boundary  value  problems  may  be  obtained  using  finite  element  approximations. 
Finite  element  methods  will  be  discussed  in  the  following  two  sections.  For 
general  references  to  finite  element  techniques,  the  reader  is  directed  to 
Huebner  (1975)  and  Zienkiewicz  (1971). 

Let  us  now  consider  how  the  problem  of  classical  plane  elastostatics  can 
be  generalized  to  encompass  adaptative  elastic  materials.  We  will  first  note 
that  the  requirements  of  symmetry  about  the  X3  axis  will  require  the  trabecular 
adaptation  tensor  M to  have  the  form 
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(144) 


i.  e.  there  exists  a trabecular  plate  normal  to  the  x-^-plane.  In  order  to 
two-dimensionalize  the  governing  equations  (122-126)  one  must  first  two- 
dimensionalize  the  tensor  M in  a manner  analogous  to  the  reduction  of  the 
stiffness  tensor  (139)  . This  will  require  that  there  exists  a function  H , 
mapping  the  set  of  all  three-dimensional  second-rank  tensors  in  the  form  of 
(144)  onto  the  set  of  all  two-dimensional  symmetric  second-rank  tensors  M , i. 

H ( M ) = M (145) 


where  this  mapping  is  subject  to  the  constraint  that 
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for  all  M 

33 


(146) 


Having  defined  such  a mapping,  the  governing  equations  (122-126)  take  the 
form 

• A 

“ab  = *ab  1 E > -P  > 

[w*>  Ecd]  ,b  = ° 

and  (147) 

2 E12,  12  = Ell,  22  + E22, 11 

and  are  subject  to  the  initial  conditions 

( x , t).  = M°  (x)  (148) 

t=  0 
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and  the  boundary  conditions 
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(x.t)  = CabcdtM(2't>^  E„H(2’t)n 


cd' 


(149) 
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while  the  stiffness  and  remodeling  functions  are  subject  to  the  constraints. 


C , ,[H(M)]  = C , .(M)  - Cab33(-)  Ccd33(~* 

abed  - - abed  - E^Tm) 
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respectively,  where 
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Let  us  now  consider  the  physical  significance  of  these  two-dimensional 
approximations.  The  existence  of  the  mapping  H(M)  implies  that  for  any 
three-dimensional  trabecular  structure  whose  trabecular  adaptation  tensor  can 
be  described  by  (144),  there  must  exist  an  equivalent  "reduced"  trabecular 
structure  in  which  the  trabeculae  normal  to  the  x3  axis  have  been  removed, 
i.  e.  , one  for  which  the  trabecular  adaptation  tensor  has  the  form 


“u 

M12 

0 

S12 

^22 

0 

(153) 

0 

0 

0 

The  constraint  on  the  stiffness  function  (150)  implies  that,  for  the  "reduced" 
trabecular  structure  • The  remaining  two  sets  of  trabeculae  which  run  along  the 
x 3 axis  have  been  strengthened  appropriately  so  that  the  stress-strain 
relationship  of  the  original  three-dimensional  trabecular  structure 
is  equivalent  to  the  stress -strain  relationship  of  plane  stress. 

Likewise,  the  constraint  on  the  remodeling  function  (151)  implies 
that  the  "reduced"  trabecular  structure  remodels  in  a manner  exactly 
equivalent  to  that  of  the  original  th  ree -di  me  nsion  trabecular  structure, 
when  the  latter  is  subjected  to  plane  stress 
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At  this  point  we  should  make  the  obvious  comment  that  the  assumptions 
necessary  in  these  two-dimensional  approximations  of  the  governing  equations 
will  not  necessarily  be  satisfied  exactly  in  the  case  of  the  living  bone  that  we 
wish  to  consider.  We  might  expect  that  the  two-dimensional  governing  equa- 
tions will  describe  behavior  different  from  that  occurring  in  the  three- 
dimensional  living  bone,  whenever  and  wherever  one  of  the  simplifying 
assumptions  discussed  above  is  violated.  For  example,  the  violation  of  the 
constraints  (150)  and  (151)  should  become  most  manifest  in  a region  in  which 
there  exist  large  variations  in  the  value  of  M33.  In  the  case  of  the  human 
femur,  which  will  be  discussed  in  subsection  III.  6,  the  M33  component  will 
vanish  in  an  anatomical  region  called  "Ward's  Triangle".  It  is  in  the  region 
of  Ward's  Triangle  that  the  two-dimensional  governing  equations  developed 
here  will  prove  least  satisfactory  in  describing  the  trabecular  adaptation  of 
the  human  femur. 

III.  3 NUMERICAL  TECHNIQUES 

a.  General  Plan  of  Solution 

We  have  already  stated  that  adaptive  elasticity  may  be  regarded  as  a 
generalization  of  classical  elasticity  in  which  the  elastic  modulii  are  considered 
to  be  functions  of  the  strain-history.  This  property  has  an  important 
application  in  the  case  of  an  elastostatic  boundary  value  problem.  Given  that 
the  boundary  conditions  are  not  a function  of  time,  the  solution  of  (147)2  and 
(147)3  for  the  strain  field  E at  the  time  t=0  can  be  readily  obtained  by  setting 
the  field  M equal  to  the  values  given  by  (148)  and  then  considering  the  solution 
of  equations  (147)2  and  (147)3  for  the  field  E as  equivalent  to  the  solution  of 
the  classical  plane  elastostatic  boundary  value  problem  (141).  Having  obtained 
the  strain  field  E , the  change  in  the  trabecular  configuration  in  some  time 
interval  A t can  be  calculated  by  numerical  integration  in  time  of  the  governing 
equation  containing  remodeling  function  (147)^.  The  value  of  the  trabecular 
adaptation  at  the  time  t+At  can  then  be  used  to  solve  (144)2  and  (144)^  for  the 
strain  field  E at  time  t+At  and  this  process  can  be  continued  indefinitely.  In 
this  way  one  can  reduce  the  problem  of  two  dimensional  adaptive  elastostatics 
into  two  simpler  problems:  (1)  solve  equations  (147)2  and  (147)3  to  determine 
the  strain  field  as  a function  of  space,  and  (2)  use  Picard  integration  of  (147  )j 
to  determine  the  trabecular  adaptation  field  M as  a function  of  time. 

b.  Solution  for  the  Strain  Field 

The  solutions  of  equations  (147)2  and  (147  >3  can  be  effected  by  the  method 
of  finite  element  approximations.  In  its  simplest  form  the  finite  element 
method  resolves  to  the  following  technique:  We  divide  the  body  into  a system 
of  v triangular  elements  at  whose  vertices  reside  a total  of  N nodes.  Next 
wc  approximate  the  governing  equations  by  assuming  that  the  strain  and  the 
physical  properties  of  the  material  within  each  of  these  V elements  are 
uniform.  Since  the  strain  and  stiffness  of  each  element  are  uniform,  it  follows 
that  the  stress  within  each  element  is  uniform.  It  is  possible  to  resolve  the 


internal  stress  of  the  *-th  element  as  three  two-dimensional  forces  (f*X,  ffy), 

(fj  , f j 7 ) and  (f^.  , f k 7 ) acting  at  the  node  vertices  i,  j,  and  k,  respectively. 

In  a similar  manner,  it  is  possible  to  resolve  the  strain  on  the  *-th  element 

x y x y x y 

in  terms  of  the  three  vectors  (u^  , u[  ) , (uj  , u:  ) and  (uk  , uk  ) , which  describe 
the  two-dimensional  displacement  of  the  three  vertices  of  the  <c-th  triangular 
element.  Using  the  stress-strain  relations,  it  is  then  possible  to  write  the  force- 
displacement  relation  for  the  *-th  element  in  the  form 
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The  body  is  nonaccelerating  and  each  of  its  component  nodes  is  nonaccelerating. 
Consequently,  the  total  effective  force  applied  to  each  node  must  be  zero. 

The  total  force  of  the  i-th  node  equals  the  total  of  the  equivalent  forces 
arising  from  the  internal  stress  in  each  of  the  elements  contingent  on  the  i-th 
node,  (f^X,  f^y),  plus  the  forces  from  the  boundary  traction  applied  to  the 
i-th  node.  ( \Ve  recall  that  we  have  assumed  that  the  body  force  is  zero.  ) If 
we  write  the  boundary  traction  force  on  the  i-th  node  with  the  array  -(T,  fy), 
then  nonaccleration  of  the  i-th  node  implies  that  1 1 

f x = 2 f 


If  one  then  uses  linear  superposition  of  equation  (154)  with  equation  (15  5),  one 
can  then  write  the  boundary  forces  in  terms  of  the  displacements  by  the 
linear  relation 
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where 
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However,  the  boundary  traction  forces  (f^  , fj7 , ...  f^,  ) can  be  resolved 

in  terms  of  the  boundary  tractions  described  in  equation  (149).  Consequently, 
using  finite  element  approximations  one  can  reduce  the  solution  of  equations 
(147)2  and  (147  >3  for  E into  the  solutions  of  the  linear  equation  (156)  for  the 
displacements  (upvj ),  (u^.v^), (UN’VN^*  Having  determined  the  dis- 

placements, the  strain  field  E can  be  calculated  directly. 


Before  we  consider  the  solution  of  (156)  for  the  displacement,  two 

y y X y 

comments  are  in  order  with  respect  to  the  arrays  (f^  , f£  f^.f^)  and 

(uf,uf . . ,uX,  Uy).  First  we  will  observe  that  the  boundary  forces 

depicted  by  (f^, fy, ,f^,f^)  are  not  linearly  independent.  Since  the 

body  is  in  static  equilibrium,  the  sum  of  all  external  (boundary)  forces  must 
equal  zero,  thus 
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Also  the  sum  of  the  moments  about  some  arbitrary  point  (x0,yQ)  must  be 
zero,  thus 

N 

0 = L [ (y.  - yo)  £ - -(x.  - xQ)  f[  ] (159) 

i = 1 

where  the  position  of  the  i-th  node  is  given  by  the  array  (x^,  y^).  Consequently, 
in  programming  any  solution  of  (156),  one  must  install  error-detecting  sub- 
routines that  will  verify  the  validity  of  constraints  (158)  and  (159).  Otherwise 
boundary  data  violating  these  constraints  can  be  inadvertently  entered,  causing 
the  program  to  yield  spurious  results. 

The  second  observation  is  that  the  solution  of  the  boundary  value  problem 
described  by  equations  (155)  through  (159)  will  not  be  unique.  It  can  be  demon- 
strated that  the  solution  of  equation  (157)  will  be  indeterminate  within  a uniform 
translation  and  a rigid  rotation.  In  order  to  render  the  solution  of  (157)  unique, 
we  must  make  the  following  two  assumptions:  First,  let  us  constrain  transla- 
tion by  assuming  that  one  of  the  internal  nodes  (i.  e.  j a node  which  is  not  on 
the  boundary)  has  a zero  displacement  and  the  ordering  of  the  nodes  is  such 
that  this  stationary  node  has  the  index  of  N.  This  will  mean  that 


and  (160) 

uy  = 0 

N 


Secondly,  we  will  constrain  rotation  by  assuming  that  one  other  internal  node 
can  move  only  in  the  x-direction  and  that  this  second  constrained  node  has  the 
index  of  N-l.  This  will  mean 


It  is  assumed  that  there  are  no  body  forces  and  that  the  nodes  of  indi  ces 
N-l  and  N are  not  on  the  boundary.  This  will  mean  that 


(162) 


When  these  constraints  and  the  constraints  (160)  are  substituted  into  (156)  we  will 
have  the  2N  equations  in  2N-3  variable,  which  can  be  expressed  by  the  matrix 
relation 
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where  A=  N-2  and  m=N-l  . Of  course  the  array  (ui  , u,  . 

’ 1 1 ' K l ' 

is  over-determined  and  it  is  notdificult  to  demonstrate  that  the  last  three 
equations  of  (163)  can  be  obtained  from  the  first  2N-3  using  Newton's  third 
law  and  the  conservation  of  momentum  equations  (158)  and  (159).  Consequently, 
the  2N-3  equations  given  by  the  matrix  equation 


x 

cY 


1 


f* 

I 


A** 

11 

AXy 

A11 

Axx 

1A 

AXy 

1A 

A 

A11 

A77 

A11 

AXy 

aia 

Ayy 

1A 

A 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

AXy 

aai 

• 

• 

• 

XX 

aaa 

AXy 

AA 

A 

Ayy 

Aii 

• 

• 

• 

AXy 

aaa 

Ayy 

AA 

A' 

XX 

‘lm 

lm 


xx 


xy 


XX 

L Am 
xy 
Am 
xx 


m A in  A mml 


u 


u 


u 


A 

•y, 

X 

u 

m 


(164) 


xy 

ml 


k 


will  uniquely  determine  the  array  (uj,uj,  ....  um  ).  This  equation  can 
be  solved  by  conventional  numerical  methods  such  as  Gauss-Seidel  elimination 

7 XV  X 

and  once  the  array  (u^.u^,  . . . , um)  is  determined,  appending  equations 
(160)  and  (161)  gives  a complete,  unique  solution  for  the  displacement  array 

(uj.Uj,.  . , , , Uj^  ) . 

c . Solution  for  the  Trabecular  Adaptation  Field 

Assume  that  we  are  given  the  trabecular  adaptation  field  M(x,  t)  for 
t^t^  and  the  strain  field  E(x,  t)  during  the  time  interval  Integra- 

tion of  the  conservation  of  mass  equation  (147)  gives  us  the  integral  equation 


M(x,t  .)  - M(x,t  ) = f R[E(x,  t),  M(x,  t ) ] dt' 
n+1  ~~n  .A  . ~ ~ 


This  integral  equation  suggests  that  if  the  boundary  tractions  are  fixed  and  the 
strain  field  E(x,  t)  can  be  considered  time -invariant  during  the  interval 
(tR+  j>tn)  > then  one  might  be  able  to  approximate  the  trabecular  adaptation 
function  at  t=tn+  \ by  the  expression 


~(~’tn+l)  ' ^(^’tn)  = SC  E(x,tn),  M(X)tn)]  Atn  (166) 

where  At  = t - t 
n n+  1 n 


and  then  use  the  techniques  of  the  previous  subsubsection  to  calculate  the  strain 
field  at  the  time  t»tn^_  ^ , Recursion  of  these  two  approximation  techniques  can 
then  be  used  to  calculate  the  strain  and  the  trabecular  adaptation  fields  at 
any  arbitrary  later  time  increment. 

In  addition,  these  numerical  methods  also  allow  one  to  consider  the 
behavior  of  an  adaptive  elastic  material  under  even  more  general  conditions. 
For  instance,  from  the  discussion  in  subsubsection  I.  4.  d,  one  might  conclude 
that  it  might  be  useful  to  study  the  properties  of  an  adaptive  elastic  material  in 
the  case  when  the  boundary  tractions  are  cyclic  and  the  function  s(x,  t)  as 
used  in  equation  (149)  is  periodic,  i.  e. 

s(x,t)  - s(x,  t+r ) for  all  xe5B  (167) 
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If  numerical  methods  are  employed,  one  of  the  simplest  periodic 
functions  which  can  be  studied  will  be  one  in  which  the  load  cycle  of  duration  r 
is  divided  into  two  hemicycles  of  equal  duration  during  which  the  loading 
assumes  the  values  of  two  time -invariant  boundary  tractions  sl^)  and 
s2(x)  , i.  e. 


s(x,  t) 

where 


s*(x)  u(t-t  ) + 

co  O 

u(t)  = s (-1) 

m=  U 


2 

s 

-m 


(x)  u(t-t  -z—  ) 
o 2 

u [ 2(t-mr)  ] 


(168) 

(169) 


and  u(t)  is  the  Heavyside  unit  step  function.  If  we  assume  that  the  inertial 
terms  in  the  conservation  of  mass  equation  can  still  be  ignored,  the  quasi-static 

ection  resolved  into  the  solution  of  the  two  strain 
, corresponding  to  the  strain  in  the  adaptive 
rabecular  adaptation  field  M(x,  f^)  and  the  boundary 
traction  fields  sL(x)  and  s4(x)  , respectively.  One  might  then  argue  that 
integration  of  the  conservation  of  mass  equation  allows  us  to  approximate  the 
trabecular  adaptation  at  the  time  t-tn+  i by  the  expression 


problem  of  the  previous  subs 
fields  .g  (x,  tn)  an<^  E (x,  t^) 
elastic  material,  given  the  t* 


M(x,  t L J - M(x,  t ) = 1/2 

~ ~ n+  1 


JrEe^x,  tn) , 


M(x,  t ) ] + R(E  (x,  t ),  M(x 


x,  t )|  At 
n I n 


(170) 


In  closing,  we  might  comment  as  to  the  final  form  the  trabecular  adapta- 
tion function  field  M(x,  t)  will  take  as  t-*00  . When  bone  is  subjected  to  static 
or  cyclic  loading  conditions,  one  of  two  things  can  be  observed  to  happen 
physiologically:  (1)  the  bone  will  remain  physiologically  stable  and  asymptoti- 
cally enter  a condition  of  homeostasis,  or  (2)  the  bone  will  become  physiolog- 
ically unstable  and  enter  a pathological  condition  (such  as  overstrain  necrosis) 
leading  to  structural  failure.  In  terms  of  the  model,  structural  failure  of 
bone  corresponds  to  the  strain  or  trabecular  adaptation  fields  attaining  values 
outside  of  the  normal,  nonpathological,  physiological  range.  Such  conditions 
of  instability  can  be  described  quite  rigorously  in  mathematical  terms  using 
the  concepts  of  mathematical  stability  introduced  by  Lyapunov.  The  problem 
of  stability  and  concomitant  structural  failure  was  the  subject  of  considerable 
discussion  in  the  author's  dissertation. 
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The  condition  of  homeostasis  corresponds  to  the  conditions  which  can 
be  expected  of  bone  with  a normal  remodeling  function  subject  to  a normal 
loading  configuration.  Under  conditions  of  homeostasis,  the  trabecular  adapta- 
tion field  M will  have  attained  values  such  that  the  ambient  loading  will 
cause  exactly  the  right  amount  of  strain  so  that  average  value  of  the  remodeling 
function  R(E,  M)  over  the  loading  cycle  is  zero  for  every  point  in  the  body. 

This  trabecular  configuration  can  be  termed  the  homeostatic  trabecular  con- 
figuration. For  example,  the  trabecular  configuration  illustrated  in  Figure  4 
corresponds  to  the  homeostatic  trabecular  configuration  for  a human  femur 
sxibject  to  normal  loading  on  a human  femur.  Determination  of  the  homeo- 
static trabecular  configuration  for  a given  remodeling  function  and  load  con- 
figuration can  be  considered  the  single  most  important  computational  result  of 
the  computer  simulation  of  an  adaptive  elastic  material. 

III.  4 APPLICATION  OF  THE  CORTICAL  BONE  MODEL  IN  THE  EVALUATION 

OF  THE  CONSTITUTIVE  FUNCTIONS 

a.  Introduction 

As  mentioned  previously,  the  distinct  disadvantage  of  numerical  methods 
is  that  in  order  to  solve  the  governing  equations  that  simulate  an  adaptive 
elastic  material  one  must  first  have  explicit  formulations  for  the  constitutive 
functions.  Unfortunately,  quantitative  experimental  data  for  such  evaluation 
is  nonexistant  since  there  has  never,  heretofore,  existed  a coherent  mathe- 
matical model  into  which  such  data  might  be  integrated.  However,  the  author,  in 
his  dissertation,  has  successfully  described  the  remodeling  properties  of 
c ortical  bone  using  well-established  qualitative  observations  to  obtain  useful 
constraints  on  the  remodeling  function. 

In  the  following  two  subsubsections,  we  will  summarize  the  basic 
features  of  the  cortical  bone  model  as  discussed  in  the  author's  dissertation. 

After  that  we  will  apply  this  model  to  the  evaluation  of  the  constitutive  functions 
of  cancellous  bone. 

b.  Cortical  Bone:  The  Basic  Model 

In  this  subsubsection,  results  of  the  author's  dissertation  will  be  re- 
viewed and  summarized.  In  this  summary,  it  will  be  convenient  to  modify 
some  of  the  notation  and  terminology  from  that  which  appears  in  the  disserta- 
tion, not  only  for  brevity  and  simplicity  of  discussion  but  also  so  as  not  to 
conflict  with  the  notation  and  terminology  used  in  the  current  discussion. 

The  mechanical  properties  of  either  cortical  or  cancellous  bone  are 
essentially  identical  to  the  mechanical  properties  of  the  extracellular  material 
composing  the  cortical  or  cancellous  bone.  In  cancellous  bone  this  matrix 
structure  roughly  approximates  orthogonal  sets  of  trabecular  plates.  In 
cortical  bone,  this  matrix  structure  consists  basically  of  cylindrically  shaped 
osteons  whose  long  axes  run  parallel  to  the  long  axis  of  the  bone  itself.  As 
cortical  bone  "adapts"  to  its  strain  environment,  the  porosity  and  radiographic 
opaqueness  of  these  osteons  can  be  observed  to  change  (Amtmann,  1971).  In 
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his  dissertation,  (Hegedus,  1976)  the  author  simulated  the  stress -adaptation  of 
cortical  bone  by  assuming  that  this  adaptation  could  be  simulated  by  the  varia- 
tion of  a single  scalar  parameter  y analogous  to  one  of  the  zero -strain  mass 
density  parameters  introduced  in  Section  II.  Let  us  refer  to  the  parameter  y 
as  the  osteonal  adaptation  parameter  and  let  us  assume  that  its  magnitude 
has  been  normalized  so  that  its  values  must  lie  within  the  unit  interval,  i.  e. 
ye  (o,  1). 

Under  these  circumstances,  the  isothermal,  elastostatic  governing 
equations  for  cortical  bone  consist  of  a conservation  of  mass  equation, 

y = c (E  , y)  (Hi) 


a stress-strain  relation 


/\ 

Tij  " CijkX  (y  ) Ek/ 


and  a conservation  of  momentum  equation 

T =0 

ij»  J 

and  a compatibility  equation 

E..  , .+  E,  . - E..  . , - E.  . . = 0 

ij»  kjl  kX,ij  ik,  j X jX,  lk 

that  are  subject  to  initial  conditions  on  the  osteonal  adaptation 

y(x,  t)J=  y°(x) 

and  the  boundary  conditions  on  the  boundary  tractions 

T. . (x,  t)  = s n 
-U  i 3 

xi^B 


For  further  development  it  will  be  useful  to  point  out  the  existence  of 
an  auxiliary  constitutive  function,  the  compliance  function  K (y),  which  can 


be  formulated  in  lieu  of  the  stiffness  function  C„k^(y).  For ^ 
stiffness  function  y , the  value  of  the  compliance  function  K. 


jkX 


value  of  the 
(y)  will  be 


defined  as  the  fourth  rank  tensor  inverse  of  the  fourth  rank  tensor  value  of  the 
stiffness  function  (y  ) . This  will  mean  that  the  compliance  function  will 

allow  the  strain  components  to  be  expressible  as  linear  combinations  of  the 
stress  components,  i.  e. 


E.. 

ij 


V/W 


kX 


(172) 
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In  the  dissertation,  the  stress -adaptation  of  cortical  bone  was  investi- 
gated in  the  special  case  of  uniaxial  compression  of  magnitude  P(t)  applied 
to  an  initially  uniform  hollow  cylinder.  These  boundary  and  initial  conditions 
very  roughly  approximate  the  conditions  prevailing  on  a bone  such  as  the  dia- 
physis  of  the  tibia  during  normal  bipedal  activity  such  as  walking.  Given  that 
the  material  (and  therefore  the  osteonal  adaptation)  is  uniform,  i.  e. 
y(x  , 0)  = y°  , the  governing  equations  (171 ) 2 through  (171  )4  are  satisfied  by 
spatially  uniform  stress,  strain  and  volume  fraction  fields  which  have  the  form 


and 


-P">  6j3 


Kljkl(r)  Pit) 


y (x,  t)  = "y  (t) 


where  y (t)  is  the  solution  to  the  equation 


y - g (P(t)  , y ) 


(173) 


(174) 


where  y (t)  is  subject  to  the  initial  conditions 


y (t)|  = y° 

t=  0 


(175) 


and  where  the  stress  remodeling  function  g(P,  y)  is  defined  by  the  relation 

A ^ 

g(P,  y)  =,  c(  -K.j33(y)P,  y)for  all  P and  y (176) 

c.  Cortical  Bone:  Phenomenological  Constraints 

At  this  point  in  the  development  we  are  left  in  a situation  somewhat 
similar  to  that  encountered  in  the  previous  subsection  in  the  analysis  of 
cancellous  bone:  a constitutive  function  of  two  variables,  for  which  there 
exists  no  quantitative  data.  However,  there  do  exist  a number  of  well-documented 
qualitative  descriptions  of  the  phenomena  of  bone  remodeling  from  which  one 
can  derive  constraints  for  the  stres  s - remodeling  function,  g(P,  y)  . These 
phenomenological  constraints  can  be  summarized  as  follows: 

(1)  There  Exists  a Compressive  Stress  for  Which  the  Rate  of 
Remodeling  is  Zero 

Matrix  material  is  constantly  being  deposited  and  resorbed, 
but,  under  normal  loading  conditions,  the  net  osteonal  adaptation 
does  not  change.  The  normal  ambient  loading  of  cortical  bone 
is  compressive,  thus  we  will  postulate  the  existence  of  some 
ambient  compressive  stress  Po>0  » such  that  with  a given  osteonal 
adaptation  yQ  , the  rate  of  remodeling  is  zero,  i.e. 

° = g(P0’V  (177) 
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(2)  A Negative  Slope  with  Respect  to  yQ  at  the  Point  of 
Zero  Remodeling. 

The  most  important  feature  of  the  remodeling  of  living 
bone  is  its  ability  to  maintain  homeostasis  stably.  One  can  assume 
that  the  application  of  the  compressive  stress  PQ  will  cause  the 
adaptive  elastic  material  to  remain  stable  at  the  osteonal  adapta- 
tion yQ  . Thus  if  y were  to  increase  slightly  from  yQ,  g(PQ,  yQ) 
must  become  negative  in  order  to  force  the  osteonal  adaptation  to 
decrease.  Likewise,  if  yQ  were  to  increase  slightly  from  yQ  , 
g(PQ,  yQ)must  become  negative.  Mathematically,  this  condition 
of  stability  requires  that 

|^g(Po,y)|<  0 (178) 


y=  y 

We  will  digress  for  a moment  from  our  discussion  of  the  phenomological 
restrictions  of  the  modeling  function  in  order  to  introduce  some  important 
definitions.  In  general,  we  will  define  a combination  (P^.  >h)  ax^a^  stress 
P^  and  osteonal  adaptation  y^  such  that 

g(ph’V  = 0 

-^-g(Ph,y")  < 0 (179) 

/=  yh 

as  a point  of  homeostasis.  The  particular  point  of  homeostasis  consisting  of 
the  volume  fraction  yQ  and  the  axial  loading  PQ  will  be  called  the  homeostatic 
datum. 

Under  normal  loading  conditions  living  bone  reacts  to  an  increase  in 
loading  by  becoming  stiffer  and  less  porous.  A point  of  homeostasis  consisting 
of  a compressive  stress  P^  and  an  osteonal  adaptation  y^  will  be  said  to  be 
a normal  point  of  homeostasis  if  we  have 

8 <P  ' V > 0 (180> 

ph=P 

for  P^>  0 in  addition  to  conditions  (177)  and  (178).  For  instance,  given  that 
y ~ Vo  ' we  might  expect  that  if  the  stress  were  to  increase  slightly  fio  m 
that  of  the  homeostatic  datum,  matrix  accretion  would  occur  and  the  net  rate 
of  reaction  would  become  positive.  Similarly,  if  this  stress  were  decreased 
slightly  below  that  of  the  homeostatic  datum,  matrix  material  would  be 
resorbed  and  the  rate  of  reaction  would  be  negative. 


On  the  other  hand,  in  certain  pathological  conditions  such  as  overstrain 
necrosis,  bone  reacts  to  an  increased  load  by  becoming  more  porous,  i.  e. 

mt  8 <p  ■ V <0  ,181> 

p=ph 

for  < 0 , in  this  case,  the  material  will  be  said  to  be  at  a paradoxical 
point  of  homeostasis.  Finally,  we  should  note  that  in  the  case  of  tensile 
stress,  (P  < 0 ) the  definitions  (180)  and  (181)  for  normal  and  paradoxical  points 
of  homeostasis  will  be  reversed. 

We  will  now  return  to  our  discussion  of  the  phenomenological  constraints 
on  the  remodeling  function.  Ln  terms  of  the  definitions  given  above,  the  stress- 
r emodeling  function  g(P,  y)  has  the  following  additional  property: 

(3)  The  Homeostatic  Datum  is  a Normal  Point  of  Homeostasis. 

This  means 


g(P,  y)  | > 0 (182) 

P=  P 

o 

(4)  Disuse  Resorption  Occurring  at  y = 0 

It  is  observed  clinically  that  if  bone  of  average  porosity  is 
subjected  to  zero  loading  conditions,  bone  matrix  will  be 
resorbed.  This  means  that  unstrained  matrix  material  under- 
goes a type  of  "disuse  atrophy"  and  thus 

g (0,  y)  < 0 (183) 

(5)  The  Maximum  Rate  of  Remodeling  is  Small 

Bone  remodeling  is  a very  slow  process.  Bone  labeling 
studies  by  Frost  (1969)  indicate  typical  remodeling  rates  on 
the  order  of  months.  Even  under  extreme  conditions  such  as 
weightlessness,  two  weeks  are  required  to  effect  changes  in 
skeletal  mass  on  the  order  of  10%.  We  can  thus  conclude  that 
the  values  of  the  material  parameters  which  determine  the 
remodeling  function  have  values  such  that  the  quantity 

1 

g(P,y)  (184) 

7 

will  be  of  the  order  of  several  weeks,  i.  e.  10  seconds. 
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d.  Cancellous  Bone:  Approximation  as  a Pinned  Triangular  Structure 

Anatomically,  there  does  not  exist  any  sharp  lineof  demarcation  between 
a region  of  cortical  bone  and  a region  of  cancellous  bone.  Furthermore, 
since  the  histological  distinction  between  cortical  bone  and  cancellous  bone  is 
more  one  of  the  structural  arrangement  of  the  lamellar  bone  composing  the 
matrix,  one  can  assume  that  the  remodeling  properties,  as  indicated  by  the 
stiffness  and  remodeling  functions  of  cortical  cancellous  bone,  are  comparable 
if  not  identical.  In  fact,  cancellous  bone  might  be  thought  of  as  a structural 
rearrangement  of  cortical  bone  into  trabecular  patterns  which  optimize  the 
load-carrying  properties. 

It  was  already  stated  in  the  previous  subsection  that  the  field  equations 
(171)  describing  cancellous  bone  can  be  approximated  by  assuming  material 
uniformity  over  a finite  triangular  region.  Let  us  now  assume  that  these 
uniform  properties  can  be  simulated  by  a pinned  triangular  structure  extending 
along  the  boundaries  of  each  element  and  each  one  of  these  structural  members 
will  be  termed  a trabecular  link,  or  more  simply,  a link.  Let  u?  further 
assume  that  each  of  these  trabecular  links  is  mechanically  equiv^.  .ent  to  a 
small  uniform  portion  of  cortical  bone  whose  properties  have  been  described 
in  the  author's  dissertation  and  assume  that  the  trabecular  links  composing 
this  structure  are  of  uniform  length  and  thus  form  an  equilateral  triangular 
structure.  We  recall  that  in  the  cortical  bone  model,  the  index  of  osteonal 
adaptation  was  a single  scalar  field  y proportional  to  the  zero  strain  mass 
density  of  cortical  bone.  Assuming  that  cancellous  bone  is  composed  of  an 
equilateral  triangular  structure  of  trabecular  links  whose  properties  are 
identical  to  cortical  bone,  is  equivalent  to  assuming  that  the  osteonal  adapta- 
tion of  cancellous  bone  can  be  described  by  the  three  scalar  fields  y^,  y , 
and  y ^ , representing  the  trabecular  adaptation  of  the  trabecular  links  along 
the  trabecular  basis  vectors  w = (1,  0)  , w^  =1/2  (l,\/3  ) and  w^  =1/2  (-l,v/3), 
respectively.  Using  equation  (112)  one  can  then  express  the  effective  trabecular 
adaptation  tensor  in  terms  of  these  three  scalar  fields  by  the  relation 
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In  this  subsection  we  will  dispense  with  the  bar  to  indicate  two-dimensionaliza- 
tion  of  the  parameters  T , E , and  M. 
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The  next  step  will  be  to  consider  the  constitutive  functions.  To  facilitate 
this  discussion,  let  us  consider  two  new  physical  parameters  that  can  be  associated 
with  each  trabecular  link.  The  first  of  these  new  parameters  will  be  termed 
the  trabecular  strain,  ett,  and  will  indicate  the  material  strain  along  the 
trabecular  links  parallel  to  the  trabecular  basis  vector  w®  . In  terms  of  the 
displacement,  the  trabecular  strain  on  a trabecular  link  connecting  the 
i-th  and  j-th  nodes  of  an  equilateral  triangular  trabecular  structure  will 
have  the  form 


e 


a 


1_ 

l 


u^  ) + 


(186) 


where  l is  the  uniform  length  of  the  trabecular  links  and  the  unstrained 
position  of  the  j-th  node  equals  the  unstrained  position  of  the  i-th  node 

fy 

plus  the  vector  . Usin^;  the  kinematic  definitions  of  continuum  mechanics, 

it  is  possible  to  express  e in  terms  of  the  strain  E by  the  relation 
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The  second  new  parameter  will  be  termed  the  trabecular  stress  0 and 
will  indicate  the  tension  (or  compression)  in  the  trabecular  links  parallel  to 
the  trabecular  basis  vector  w®  . In  terms  of  the  two-dimensional  stress  T 
(whose  dimensions  are  force  per  unit  length)  the  trabecular  stress  can  be 
expressed  by  the  relation 

cy  of  fy 

0 = £w  T w (188) 


e.  Cancellous  Bone:  The  Constitutive  Functions 


We  are  now  in  a position  to  discuss  the  constitutive  functions  themselves. 
For  simplicity,  we  will  assume  that  the  stiffness  of  any  one  of  these  trabecular 
links  is  directly  proportional  to  the  osteonal  adaptation  along  the  axis  of  that 
particular  trabecular  link.  We  might  note  that  the  assumption  that  the  stiffness 
of  the  trabecular  link  is  directly  proportional  to  the  osteonal  adaptation  y® 
is  equivalent  to  assuming  that  the  material  composing  each  trabecular  link 
is  uniform  and  the  osteonal  adaptation  is  directly  proportional  to  the  effective 
thickness  of  each  trabecular  link.  This  implies  that 
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where  > is  a material  constant  with  the  dimensions  of  force  per  unit  length. 

It  is  then  possible  to  demonstrate  that  the  force-displacement  relationship  for 
the  trabecular  link  parallel  to  w®  connecting  the  i-th  and  j-th  nodes  will 
have  the  form 
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Since  the  structural  members  are  pinned,  the  force  exerted  by  each  trabecular 
link  will  depend  on  only  the  displacement  of  the  two  end  points  of  that  trabecular 
link.  Consequently,  the  force -displacement  matrix  for  the  structure  as  a 
whole  can  be  obtained  by  summing  the  expression  (190)  over  each  trabecular 
link  A in  the  structure. 

Let  us  now  consider  the  remodeling  function:  Let  us  suppose  that  the 
remodeling  occurring  in  any  given  trabecular  link  depends  only  on  the  physical 
c onditions  within  that  particular  link.  This  will  mean  that  the  rate  of  reaction 
of  a trabecular  link  will  depend  on  only  the  trabecular  strain  e a and  the 
trabecular  adaptation  7U  . Consequently,  the  conservation  of  mass  equation 
will  take  the  form 


(191) 


where  the  link  remodelinp^function  ?a.(r®y**)  can  be  expressed  in  terms  of 
the  remodeling  function  R(E,  M)  by  the  expression 
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Let  us  further  suppose  that  in  view  of  the  quadratic  dependence  of  the  remodeling 
function  discussed  in  section  II  , the  link  remodeling  function  is  a 
quadratic  function  of  both  the  trabecular  strain  and  the  trabecular  adaptation, 
i.  e.  , 

"“(«•>)"  a0+  al*  + a2y+  V2+  Vy+  a5yZ  (194) 


where  aQ  , a,  , a£  > a 3 , a4  , and  a$  are  material  constants.  We  might  also 
note  that  the  link  remodeling  function  used  in  equation  (194)  can  also  be  recast 
in  terms  of  the  link  stress  0a.  This  allows  us  to  recast  (191)  into  the  form 
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where 


• a a,„a  y. 

y = g (o  , yT ) 


g(a,y)=  aQ+  — + a2y  + a3 


(19  5) 


+ a5y  (196) 


We  can  now  assume  that  the  stress -remodeling  function  g(o,7)  is  exactly 

analogous  to  the  stress  - remodeling  function  defined  by  equation  (176),  which 

describes  the  remodeling  properties  of  cortical  bone,  except  for  the  sign  of  its  first 
argument. 

Let  us  now  consider  the  restrictions  of  g(o,7)  that  are  imposed  by  the 
phenomenological  constraints.  First  of  all,  let  us  assume  that  the  force  and 
displacement  variables  have  been  normalized  so  that  we  can  set  the  material 
constant  K equal  to  unity.  Next  let  us  assume  that  the  compressive  homeo- 
static datum,  introduced  previously  has  the  values  of  a=-land  7=1  and  hence. 


g (-1,D=  o 


(197) 


Next  we  will  observe  that  the  stability  of  the  compressive  point  of  homeostasis 
will  require  that 
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In  addition,  we  can  also  assume  that  the  compressive  point  of  homeostasis  is 
a normal  point  of  homeostasis,  hence 

1.  g (CT  , l)i  <0  (199) 

I 

a=  -1 


The  particular  nature  of  cancellous  bone  implies  further  constraint. 

It  has  long  been  accepted  that  the  trabeculae  of  cancellous  bone  are  capable 
of  maintaining  homeostasis  under  tensile  loads.  In  fact,  in  the  human  femur, 
whose  trabecular  patterns  are  illustrated  in  Figure  4,  two  of  the  five  major 
trabecular  groups  carry  tensile  loads  and  in  fact  are  referred  to  as  the 
Principal  and  Secondary  Tensile  Groups.  The  arguments  in  favor  of  tensile 
conditions  of  homeostasis  also  have  indirect  evidence.  It  is  well  known  that 
moderate  cyclic  loading  will  cause  greater  trabecular  development  than  static 
loading  of  the  same  average  magnitude.  In  his  dissertation,  the  author 
demonstrated  that  this  phenomenon  would  occur  if  and  only  if  there  existed 
both  compressive  and  tensile  points  of  homeostasis  and  the  coefficients  of  the 
remodeling  function  were  such  that 


90 


(200) 


> 


g(cr  , y)|>  o 
or=  - 1 


In  closing,  it  is  important  to  cite  on  additional  constraint  on  the  remodel- 
ing function  of  cancellous  bone,  which  was  not  mentioned  in  our  earlier 
discussion  of  the  phenomenological  constraints  for  cortical  bone:  We  will 
require  that  the  remodeling  function  of  cancellous  bone  must  depend  on  the 
trabecular  adaptation  tensor  M , i.  e. 

R (E  , M)  ^ 0 . (201) 


The  proof  of  this  constraint  can  be  summarized  as  follows:  As  mentioned 
previously,  normal  bone,  subject  to  normal  loading  conditions  can  always  be 
observed  to  remodel  in  such  a way  as  to  asymptotically  tend  to  a condition  of 
homeostasis.  We  might  further  expect  that  such  a condition  of  homeostasis 
would  be  attainable  under  any  loading  configuration,  provided  that  loading 
configuration  were  within  the  physiological  bounds  for  that  particular  bone.  In 
particular,  a static  load  will  cause  bone  to  asymtotically  approach  conditions 
of  homeostasis  for  which  the  strain  and  the  trabecular  adaptation  approach 
time -invariant  homeostatic  values  E°  and  M°  , such  that 

R (E°  , U°)  = 0 . (202) 


For  a two-dimensional  strain  field,  E has  only  three  degrees  of 
freedom.  If  inequality  (201)  is  not  valid,  and  the  remodeling  function  I* 
depends  only  on  E , the  conditions  of  homeostasis  (202)  will  cause  the  field 
E to  be  subject  to  three  constraints.  Consequently,  if  inequality  (201)  is 
not  valid,  the  strain  field  at  homeostasis  E°  can  only  have  a finite  number 
of  values.  However,  the  material  must  also  satisfy  the  conditions  of  com- 
patability,  which  in  two  dimensions  take  the  form 


Ell,  22  + E22, 11  2 E12,  12  ‘ 


(203) 


Obviously,  the  compatibility  equation  (203)  and  the  conditions  of  homeostasis 
cannot  be  satisfied  simultaneously  if  the  boundary  conditions  arc  such  that 
the  strain  cannot  be  uniform. 

In  retrospect,  this  constraint  appears  obvious.  However,  if  one  were 
to  attempt  an  ad-hoc  computer  simulation  of  cancellous  bone  without  thoroughly 
investigating  the  mathematical  nature  of  this  simulation,  there  would  be  a 
strong  temptation  to  overlook  dependence  of  the  trabecular  adaptation  in  the 
remodeling  function  and  perhaps  a tendency  to  not  fully  appreciate  the  signif- 
icance of  the  compatibility  equation.  There  is  in  the  literature  some  indirect 
evidence  that  this  oversight  has  in  fact  occurred. 
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III.  5 THE  PROGRAMMING  PACKAGE 
a.  Introduction 


In  this  subsection  we  will  discuss  the  organization  and  application  of  a 
programming  package  that  simulates  an  adaptive  elastic  material  using  the 
numerical  techniques  introduced  in  the  three  previous  subsections.  This 
programming  package  was  designed  for  use  on  the  SCOPE  operating  system 
(a  time -shared  CDC  6600),  which  was  available  at  the  Aerospace  Medical 
Research  laboratory  (AMRL)  Initial  input  data  and  program  control  informa- 
tion are  punched  in  on  standard  80  column  cards  read  through  a Model  30  5 High- 
Speed  Card  Reader  located  at  the  AMRL  computer  terminal.  Program 
output  as  well  as  diagnostic  error  messages  are  printed  out  on  an  off-line 
I ~ charac te r line  printer  which  is  also  located  on-site.  The  program  in- 
structions are  written  in  FORTRAN  EXTENDED  VERSION  4 and  utilize  the 
standard  scientific  subroutine  package  which  is  built  into  the  SCOPE  compiler. 

In  its  output  operations,  the  programming  package  can  also  employ  an  11-inch 
CALCOMP  incremental  plotter  which  is  also  located  at  the  computer  terminal. 
The  use  of  the  CALCOMP  plotter  requires  implementation  of  auxiliary  software 
which  is  stored  in  permanent  files  on  the  SCOPE  operating  system. 

As  developed  previously,  an  adaptive  elastic  system  is  approximated 
as  an  equilateral  triangular  structure  with  N nodes.  The  main  limitation  of 
the  present  programming  package  to  simulate  this  structure  lies  in  its  ability 
to  solve  eq.  (164)  for  the  nodal  displacement.  This  solution  requires  the 
storage  of  a (2N  x 2N)  array  representing  the  force-displacement  relation. 

On  the  SCOPE  operating  system  the  maximum  total  effective  core  memory 
available  is  165K(octal)  for  daytime  runs  (800-1600  hrs.  ),  250K  for  evening 
runs  (2nd  shift  from  1600  to  2400  hrs.  ),  and  310K  for  night  runs  (3rd  shift 
from  2400  hrs.  to  0800  hrs.  ) . The  memory  capacity  of  310K  has  the  decimal 
equivalent  of  102400  or  (320)^  the  computer  could  conceivably  store  a force- 
displacement  array  for  a 160-node  structure.  However,  since  some  of  the 
memory  must  be  used  to  store  also  the  programming  instructions,  the  practi- 
cal limit  of  nodal  capacity  was  found  to  be  128  nodes  for  over-night  runs  and 
80  nodes  for  daytime  runs. 

b.  Organization 

In  order  to  reduce  the  memory  requirement  that  must  be  used  for 
programming  instructions  during  the  critical  phase  of  the  calculations  (the 
solution  of  (164))  the  calculation  sequence  has  been  broken  up  into  a series 
of  distinct  programs  that  must  be  executed  sequentially,  with  intermediate 
results  stored  on  magnetic  tapes.  This  sequence  of  programs  will  be  called 
a programming  package  and  it  will  have  the  following  four  components: 


date  nonuniform  values  of  M,  but  redesigning  the  input  format  of  INPUTR 
to  accommodate  nonuniform  values  of  M would  render  its  already  complex 
input  data  configuration  even  more  complicated.  To  accomplish  the  same 
effect,  we  can  employ: 

(4)  An  Auxilary  Program  Called  MODIFIR 

The  program  MODIFIR  allows  us  to  consider  nonuniform 
initial  values  of  M using  a relatively  simple,  fool-proof 
technique:  Initial  values  of  the  trabecular  configuration  are 
obtained  f k>  m taped  records  of  the  final,  homeostatic  trabecular 
configuration  in  a previous  calculation.  A new  value  of  the 
remodeling  function  can  be  read-in,  if  desired,  and  the  boundary 
conditions  in  the  new  calculation  are  read-in  as  scalar  multiples 
of  the  bcundary  conditions  of  the  old  calculation.  The  output  of 
the  new  calculation  is  recorded  on  logical  tape  unit  12,  using 
exactly  the  same  data  format  that  had  been  used  on  logical  tape 
unit  6 for  the  old  calculation. 

This  uniformity  of  output  between  the  output  of  MAIN 
and  MODIFIR  has  two  advantages.  First,  it  allows  the 
direct  use  of  the  output  program  RERUN  to  display  the  results 
of  the  calculations  in  MODIFIR.  Secondly,  it  allows  the  final 
trabecular  configuration  with  another  calculation  with  MODIFIR 
using  yet  another  remodeling  function  and/or  boundary  con- 
figuration. 

c.  Application:  Simulation  of  the  Stress -Adaptation  of  the  Femur 

In  this  subsection  we  will  consider  the  application  of  the  programming 
package  to  simulate  the  stres s -adaptation  of  the  human  femur  in  situ.  As 
mentioned  previously,  the  programming  package,  in  its  present  state  of 
development,  is  capable  of  calculating  the  adaptation  of  a structure  of  only 
128  nodes.  Figure  11  illustrates  the  geometric  configuration  that  was  used 
to  approximate  the  shape  of  the  human  femur.  This  configuration  was  selected 
as  the  best  compromise  between  the  conflicting  desires  to  include  as  much  as 
possible  of  the  distal  end  of  the  femur,  while  representing  as  fine  a structure 
as  possible  for  the  proximal  end  of  the  femur. 

The  first  set  of  calculations  to  be  considered  calculated  the  trabecular 
adaptation  of  an  initially  homogeneous  model  femur.  In  the  program  INPUTR 
the  remodeling  function  took  the  form 
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(1)  An  Initializing  Program  called  INPUTR 

INPUTR  reads  in  the  data  from  the  punched  cards  and 
issues  diagnostic  messages  if  it  detects  any  errors  on  the 
punched  data  input.  For  instance,  it  will  print  out  a warning 
diagnostic  message  if  the  boundary  conditions  (149)  fail  to 
satisfy  the  conservation  of  momentum  equations  (158)  and 
(159).  If  the  program  INPUTR  detects  no  explicit  error  on 
the  punch  data  input,  a predigested  form  of  the  input  data 
is  stored  on  logical  tape  unit  4. 

(2)  A Main  Iterative  Program  Called  MAIN 

The  Program  MAIN  takes  the  predigested  data  residing 
on  logical  tape  unit  4 and  performs  the  iterative  calculations 
outlined  in  Section  III.  First  it  solves  equation  (164)  for  the 
displacement  and  strain  on  each  trabecular  link  during  each 
hemicycle.  If  the  strain  on  any  one  of  the  trabecular  links 
during  either  hemicycle  exceeds  a predetermined  compressive 
or  tensile  strain  limit,  further  calculations  are  terminated 
and  an  appropriate  error  message  is  issued  by  the  line  printer. 

If  the  strain  on  each  trabecular  link  is  within  bounds,  the  con- 
servation of  mass  equation  (191)  is  utilized  to  determine  the 
appropriate  osteonal  adaptation  during  the  time  interval  At. 

The  strain  fields  are  then  calculated  for  the  next  time  increment 
and  the  process  is  repeated  for  as  many  time-increments  as 
are  called  for  by  the  input  data.  The  results  of  these  calcula- 
tions for  each  time  increment  are  stored  on  logical  tape  unit  6. 

(3)  An  Output  Program  Called  RERUN 

A third  program  takes  the  information  results  of  the 
calculations  from  MAIN  which  have  been  stored  on  logical 
tape  unit  6 and  effect  tabular  and/or  graphical  display  of 
these  calculations.  There  exist  eight  output  options,  any 
combination  of  which  can  be  selected  by  the  insertion  of 
a single  card  read  into  the  305  card  reader.  The  execution 
of  a full  optioned  graphical  output  can  be  very  time  consuming 
on  the  incremental  plotter  It  often  takes  upwards  of  20  min- 
utes to  display  a single  frame  depicting  a 128-node  structure. 

For  this  reason,  it  is  often  advisable  to  conduct  the  implemen- 
tation of  the  output  operation  in  two  parts.  First,  a cursory 
survey  to  verify  the  general  form  of  the  results  and  an  un- 
abbreviated output  only  later,  when  monopolization  of  the 
incremental  plotter  will  not  interfere  with  work  of  other 
terminal  users. 

At  the  present  state  of  development,  the  program  INPUTR  can  only  con 
sider  initial  conditions  for  the  Trabecular  Adaptation  Tensor  M which  are 
scalar  multiples  of  the  unit  tensor.  It  would  be  useful  to  be  able  to  accommo 
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where  the  material  constants  a0, iq,  a^,  eT^»  a^.  and  a^  were  read  in  as  data 
while  the  auxiliary  material  constant  Tq,  which  fix^d  time  normalization, 
was  internally  programmed  to  have  the  value  of  10  . For  the  first  set  of 
calculations  we  will  let  the  coefficients  of  the  remodeling  function  (204) 
have  the  values 


a2  _1  (205) 

and  a_  = ac  = 1 

i 5 


It  is  easily  demonstrated  that  these  coefficients  satisfy  the  constraints  discussed 
in  the  previous  subsection.  The  boundary  tractions  for  the  first  set  of  calcu- 
lations are  illustrated  in  Figure  12  and  roughly  approximate  the  boundary 
tractions  illustrated  in  Figure  5. 


The  program  calculated  that  the  model  simulating  the  femur  approached 
conditions  of  homeostasis  after  five  units  of  normalizeo  time.  The  Link 
Adaptation  of  the  model  at  homeostasis  is  illustrated  in  Figure  13.  One  of 
the  output  options  of  the  program  RERUN  allows  calculation  of  the  eigen  values 
and  eigen  vectors  of  the  trabecular  adaptation  tensor  M for  each  element  in 
the  model  structure.  The  results  of  these  eigen  value/eigen  vector  calcula- 
tions for  the  model  at  homeostasis  are  illustrated  in  Figure  14.  As  men- 
tioned previously,  the  arrangement  of  these  eigen  vector  and  eigen  values  can 
be  thought  of  as  representing  the  effective  trabecular  adaptation  of  an  equiva- 
lent orthogonal  trabecular  structure.  The  orientation  of  this  effective 
orthogonal  structure  obtained  from  the  computer  simulation  of  the  human 
femur  is  strikingly  similar  to  the  patterns  of  the  trabecular  structure  in  a 
living  human  femur  that  are  visible  in  radiographs.  Figure  15  illustrates 
this  comparison  by  copying  the  computed  compressive  trabeculae  from 
Figure  14  next  to  a copy  of  the  anatomically  observed  compressive 
trabeculae  first  shown  in  Figure  4 Likewise,  Figure  16  compares 
the  computed  tensile  trabeculae  from  Figure  14  next  to  the  anatomically 
observed  tensile  trabeculae  of  Figure  4. 


Using  the  Program  MODIFIR  it  is  also  possible  to  simulate  the 
trabecular  dissolution  that  occurs  during  disuse  and  senile  osteoporosis. 

In  Figure  17  we  have  illustrated  the  homeostatic  Link  Auaptation  that  occurs 
if  the  remodeling  function  remains  fixed  but  the  boundary  forces  are  diminished 
by  a factor  of  one-half.  These  conditions  might  simulate  the  trabecular 
dissolution  that  could  occur  when  an  astronaut  is  subjected  to  extended 
periods  of  weightlessness  or  a bedridden  convalescent  is  subjected  to  extended 
periods  of  bedrest. 
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Figure  14. 


Effective  Trabecular  Adaptation  and  Orientation  of 
Simulation  of  a ’Normal'  Femur. 
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Figure  16.  Comparison  of  the  Effective  Trabecular  Patterns  from  a Computer  Simulation  of  a 
'Normal'  Femur  and  the  Tensile  Trabeculae  of  a Living  Femur  in  Vivo. 


Clinically,  the  most  important  type  of  osteoporosis  is  senile  osteoporosis. 
Although  the  Program  MODIFIR  can  be  used  to  simulate  senile  osteoporosis, 
lack  of  time  has  so  far  prevented  this  simulation. 


Physiologically,  senile  osteoporosis  corresponds  to  a change  in  the  re- 
modeling properties  of  bone  in  such  a manner  that  a fixed  load  will  cause  the 
remodeling  function  to  become  relatively  more  negative.  In  the  author's 
dissertation  he  demonstrated  that  this  type  of  change  in  the  remodeling  function 
would  cause  bone  to  undergo  a loss  of  matrix  material  even  though  the  load 
remained  fixed. 

There  are  any  number  of  ways  that  the  remodeling  coefficients  of  (204) 
could  be  changed  so  as  to  render  the  remodeling  function  more  negative.  For 
instance,  if  the  coefficient  aQ  took  on  a negative  value,  this  might  correspond 
to  a systematic  increase  in  osteoclastic  resorption  caused  by  a homonal  change 
in  the  calcium  metabolism.  On  the  other  hand,  a decrease  in  the  coefficient  a^ 
might  correspond  to  a decreased  sensitivity  of  the  remodeling  mechanism  to 
strain.  This  decrease  in  sensitivity  might  be  caused  by  a degeneration  of  the 
piezoelectric  properties  of  collagen,  stemming  from  an  entirely  different 
hormonal  change. 

From  this,  one  might  conclude  that  what  is  presently  termed  "senile 
osteoporosis"  might  actually  be  two  (or  even  more)  distinct  pathological 
conditions,  having  different  etiologies  and  requiring  different  treatme  its. 
Differentiation  of  these  types  of  "senile  osteoporosis"  will  require  new 
diagnostic  instrumentation  that  can  objectively  measure  the  pattern  and  degree 
• of  trabecular  loss.  The  design  of  such  instrumentation  will  be  described  in 

Section  IV. 

III.  6 SUMMARY  AND  CONCLUSIONS 

In  this  section  we  have  discussed  computer  simulation  of  the  governing 
equations  developed  in  the  previous  section,  using  two  techniques:  Two- 
dimensionalization  and  finite  element  approximation.  There  are  certain 
limitations  to  these  techniques  and  to  the  computer  programs  we  have  employed. 
The  most  important  limitation  is  that  the  numerical  solution  for  the  two- 
dimensional  displacement  of  any  N node  elastic  material,  whether  it  is 
adaptive  elastic  or  classically  elastic,  will  still  require  the  solution  of  a 
linear  equation  of  order  2N  and  the  storage  of  a 2Nx2N  matrix,  representing 
the  force-displacement  relation.  It  is  obvious  that  a 128-node  representation 
of  a femur  only  crudely  approximates  the  two-dimensional  shape,  but  direct,, 
core  storage  of  its  force-displacement  matrix  will  require  approximately  10 
memory  storage  locations.  This  memory  requirement  taxes  the  storage 
capabilities  of  even  a relatively  large  computer  system  such  as  the  CDC  6600. 

There  are  a number  of  ways  that  the  effect  of  these  limitations  can  be 
reduced.  First,  it  should  be  possible  to  use  a nonuniform  mesh  in  selecting 
nodal  points.  This,  for  instance  in  the  case  of  the  simulation  of  femural 
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adaptation,  would  allow  finer  resolution  of  the  trabecular  development  in  areas 
of  special  interest,  such  as  the  proximal  end  of  the  femur,  without  wasting 
nodal  capacity  on  unirteresting  uniform  regions  such  as  the  diaphysis. 

Secondly,  it  might  be  possible  to  employ  more  economical  methods  of  matrix 
storage  which  exploit  the  symmetry  and  "bandedness"  of  the  force-displacement 
matrix.  Finally,  it  might  be  possible  to  employ  auxiliary  memory  devices, 
such  as  a high-speed  disc  memory.  The  SCOPE  OPERATING  SYSTEM  MANUAL 
alludes  to  the  existence  of  such  a feature,  called  the  ECM,  however  this 
feature  apparently  has  not  been  installed  in  the  SCOPE  system  available  at 
the  AMRL. 

Nevertheless,  albeit  limited,  the  present  programming  package  has 
demonstrated  the  feasibility  of  the  cancellous  bone  model.  It  has  predicted 
trabecular  adaptation  strikingly  similar  to  that  of  a human  femur  in  situ. 
Furthermore,  it  has  even  demonstrated  the  feasibility  of  simulating  pathological 
conditions  such  as  senile  and  disuse  osteoporosis. 


SECTION  IV 


SUGGESTIONS  FOR  FUTURE  WORK 


IV.  1 INTRODUCTION 

The  most  serious  difficulty  in  our  development  of  a theory  of  the  adaptive 
remodeling  in  bone  tissue  has  been  the  scarcity  of  quantitative  data, 
particularly  with  respect  to  the  rate  of  matrix  apposition  and  resorption  as 
a function  of  the  loading.  Consequently,  we  have  had  to  depend  exclusively 
on  qualitative  descriptions  of  bone  remodeling.  This  lack  of  data  is  understand- 
able historically.  Systematic  experiments  to  accurately  measure  the  elastic 
modulii  of  structural  materials  had  to  be  preceded  by  the  development  of  a 
mathematical  theory  of  elasticity.  The  development  of  a mathematical  theory 
of  elasticity  was,  in  turn,  motivated  by  a desire  to  systematize  experimental 
observations  of  the  behavior  of  these  structural  materials.  Likewise, 
experimental  and  clinical  observations  of  stres s -controlled  remodeling  of 
living  bone  has  motivated  a theory  of  adaptive  elasticity. 

The  quantification  of  the  stress-adaptation  properties  of  cancellous  bone 
can  be  effected  in  two  steps: 

(1)  Develop  experimental  techniques  and  instrumentation  to 

measure  the  Trabecular  Adaptation  M for  a given  sample  of  bone 

(2)  Determine  the  Constitutive  Functions  C...  .(M)  and  R(E,M). 

ijki  ~ ~ - ~ 

These  two  steps  will  be  the  subjects  of  the  following  two  subsections: 

IV.  2 MEASUREMENT  OF  TRABECULAR  ADAPTATION  TENSOR,  M 

The  concept  developed  in  this  paper  of  using  a second- rank  tensor  to 
quantify  the  trabecular  adaptation  of  cancellous  bone  is  new.  However, 
there  do  exist  in  the  literature  reports  by  workers  who  have  attempted  to 
quantify  the  trabecular  adaptation  of  cancellous  bone,  tensorial  nature  not- 
withstanding. Rockoff  f 19  68 ) gave  a semiquantitative  description  of  the  trabecular 
structure  of  a vertebral  body  using  straightforward  dissection  techniques.  Also 
recently,  Raux  et  aL  (1975)  have  given  a quantitative  description  of  the  trabecular 
structure  of  the  human  patella  using  the  techniques  of  statistical  microscopy 
described  in  Underwood  (1970).  Nevertheless,  the  dissection/microscopic 
techniques  described  by  these  authors  are  extremely  time-consuming  and 
consequently  would  render  experimental  data  reduction  prohibitively  expensive. 

We  will  now  consider  methods  of  measuring  the  trabecular  structure  of 
bone  from  radiographs.  Radiographic  measurements  are  desirable  not  only 
because  they  are  nondestructive,  but  also,  as  we  shall  see,  they  can  easily 
be  automated.  As  mentioned  previously,  the  trabecular  plates  of  cancellous 
bone  are  often  visible  on  radiographs.  They  become  particularly  clear  if 
one  of  the  trabecular  plates  runs  normal  to  the  axis  of  the  X-ray  beam.  In 
this  case,  the  trabecular  plate  normal  to  the  beam  axis  will  appear  as  only  a 
change  in  the  X-ray  density,  while  the  other  two  trabecular  plates,  which  run 
parallel  to  the  axis  of  the  beam,  would  clearly  stand  out  in  a radiograph  as 
"trabecular  lines". 
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These  patterns  of  "trabecular  lines"  , indicating  the  direction  and 
thickness  of  the  trabecular  plates,  can  be  measured  quantitatively  on  a device 
called  a polarimeter,  which  works  on  the  principle  of  Moire  interference. 

A trabecular  line  at  the  angle  0 will  be  transmitted  through  a polarimeter 
oriented  at  the  angle  6a  by  an  amount  linearly  related  to  the  quantity 
cos  2 (0  - 0a).  If  we  consider  the  transmission  factors  at  the  angles 
0l  = o°  , 02  = 120°  , and  03  = 240°  as  proportional  to  the  trabecular  adapta- 
tion parameters  y1  , y2  , and  y3  , respectively,  then  employing  eq.  (185) 
determines  the  value  of  M for  each  point  of  the  radiograph. 

Since  the  polarimeter  is  coming  into  more  frequent  use  at  a number  of 
research  institutions,  it  is  important  to  point  out  a possible  error  in  the 
experimental  protocol  of  its  use:  Some  researchers  have  attempted  to  collect 
radiographic  data  on  the  trabecular  structure  of  cancellous  bone  using  only 
two  polarimetric  readings.  In  view  of  the  tensorial  nature  of  the  trabecular 
structure  of  cancellous  bone,  data  gathered  in  this  manner  will  of  necessity 
be  incomplete.  This  error  is  exactly  analogous  to  someone  attempting  to 
completely  specify  the  strain  state  of  a plane  elastic  body  using  only  two 
strain  gauges. 

We  will  next  consider  how  the  trabecular  patterns  that  are  visible  on  a 
radiograph  and  measurable  on  a polarimeter  can  be  analyzed  mathematically. 
First  of  all,  let  us  assume  that  the  sample  of  bone  under  measurement  has 
a trabecular  structure  symmetric  about  the  axis  of  the  X-ray  beam.  If  the 
effect  of  the  trabecular  plate  normal  to  the  X-ray  beam  were  subtracted  out, 
the  resulting  radiograph  would  represent  a two-dimensional  structure  similar 
to  that  introduced  in  Section  III,  which  could  be  evaluated  by  a two- dimensional 
symmetric  tensor  Mij  . This  tensor  can  be  uniquely  resolved  as  the  sum  of 
a scalar  times  the  unit  tensor  plus  an  isochoric  tensor,  i.  e. 
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We  should  mention  two  important  properties  of  the  isochoric  tensor  f u j 
'irst,  the  eigen  vectors  of  exactly  coincide  with  those  of  the  tensor  M 


Secondly,  the  mathematical  properties  of 
multiplication  are  exactly  equivalent  to 


of  f U Vj  with  r 
those  of  the  c 


respect  to  additionaand 
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(208) 


Consequently,  if  we  consider  the  plane  of  the  radiograph  x^  - X£  to  be 
equivalent  to  the  complex  plane  x = Xj  + i X£  , the  value  of  the  trabecular 
adaptation  tensor  M can  be  uniquely  represented  by  the  real  function  of  a 
complex  variable  <f>(z ) and  a complex  function  of  a complex  variable  <f>( z ). 

It  is  not  difficult  to  demonstrate  that  the  orthogonality  of  the  trabecular 
structure  constrains  the  function  l/)(z)  to  be  an  analytic  function  of  z. 

The  two  functions  0(z)  and  0( z)  can  be  interpreted  physically  in  terms 
of  the  trabecular  structure  of  cancellous  bone.  If  we  consider  and  k2  are 
arbitrary  real  constants,  then  the  function  0(z)  maps  the  lines  z =•  kj  and 
z = ik  into  the  trabecular  lines  that  are  visible  on  the  radiograph.  We  might 
therefore  conclude  that  the  function  j/i(z)  parametrizes  the  shape  of  the  trabe- 
cular plates  visible  on  a radiograph.  Furthermore,  from  eq.  (207),  we  see 
that  the  value  of  0 is  equal  to  half  the  trace  of  M and  thus  0 also  equals 
the  average  of  the  two  eigen  values  of  M . We  might  therefore  conclude  that 
0(z)  parametrizes  the  mean  thickness  of  the  trabecular  plates  visible  on  a 
radiograph. 

We  might  note  in  passing  that  an  interesting  parallel  exists  between  our 
formulation  of  M in  terms  of  two  functions  of  a complex  variable  0 and  ij) 
and  an  analogous  formulation  of  the  stress  tensor  T used  in  complex  variable 
analysis  of  plane  elastostatics  (Sokolnikoff , 19 56 ).  In  plane  elastostatic s 
the  analog  of  the  complex  function  (z)  represents  the  shear  stress  while 
the  analog  of  the  real  function  0O  (z)  represents  the  hydrostatic  stress.  If 
we  claim  that  at  homeostasis  the  trabeculae  of  an  adaptive  elastic  material 
do  not  carry  shear  loads,  this  will  imply  that 

Arg  (til ) = Arg  (ip  ^ ) . (209) 

At  some  later  time  it  might  be  possible  to  generalize  classical  complex 
analysis  of  plane  elastostatics  so  as  to  include  the  effects  of  adaptive  elasticity. 

However,  the  formulation  of  the  trabecular  patterns  on  a radiograph  in 
terms  of  two  analytic  functions  is  more  than  simply  an  elegant  mathematical 
ftieory:  It  has  extremely  practical  implications  in  the  design  of  instrumentation. 
If  the  polarimeter  is  used  in  conjunction  with  an  image  digitalization  device, 
such  as  a flying-spot  scanner,  the  mathematical  formulation  discussed  above 
can  be  used  to  devise  computer  algorithms  with  which  data  can  be  "fit"  and 
from  which  data  "noise"  might  be  filtered.  A "smoothed"  evaluation  of  0(z) 
will  give  a continuous  quantification  of  the  mean  thickness  of  the  trabecular 
plate.  This  smoothed  data  will  not  only  allow  a much  clearer  view  of  the 
trabecular  structure,  but  more  importantly,  for  the  first  time  will  give  an 
objective  and  experimental  repeatable  quantification  of  trabecular  adaptation. 

The  design  and  development  of  an  automatic  polarimeter  will  have  clinical 
as  well  as  research  applications.  Presently,  the  most  sophisticated  method 
of  clinically  evaluating  the  degeneration  of  osteoporosis  is  a semiquantitative 
method  introduced  by  Singh  (1970)  . The  development  of  an  automatic  polari- 
meter would  give  the  clinician  a more  accurat  e and  more  objective  criterion 
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for  diagnosing  osteoporosis.  Furthermore,  we  have  already  pointed  out  that 
the  mathematical  analysis  suggests  that  there  are  more  than  one  type  of 
osteoporosis.  Tnese  different  types  of  osteoporosis  can  be  expected  to 
induce  quantitatively  different  osteoporotic  trabecular  patterns.  Therefore 
the  development  of  an  automatic  polarimeter  would  be  of  great  clinical  importance 
in  the  diagnosis  and  treatment  of  osteoporosis. 

IV.  3 MEASUREMENT  OF  THE  CONSTITUTIVE  FUNCTIONS 

There  are  two  constitutive  functions  to  be  determined:  the  Stiffness 
Function  and  the  Remodeling  Function.  Let  us  first  consider  the  Stiffness 
Function.  At  first  glance,  the  evaluation  of  the  Stiffness  Function  C-^^M) 
seems  formidable  indeed.  The  stiffness  tensor  has  21  independent 
parameters  that  depend  on  6 independent  variables.  However,  the  symmetry  and 
objectivity  relationships  discussed  in  Section  II  yield  the  formula  (115),  which 
reduces  the  formulation  of  the  stiffness  of  cancellous  bone  (as  a whole  tissue) 
down  to  the  considerably  simpler  problem  of  determining  the  stiffness  function 
of  a single  trabecular  plate.  There  are  a number  of  researchers,  among  them 
Townsend  and  Rose  (1975)  , who  are  presently  investigating  the  problem  of 
the  elastic  properties  of  a single  trabecular  plate. 

A Let  us  next  consider  the  problem  of  determining  the  Remodeling  Function 
Rjj(E,  M).  Again  the  objectivity  and  symmetry  relationships  developed  in 
Section  II  yield  a simplified  formulation  (120)  which  reduce  the  problem  of 
determining  the  remodeling  function  of  cancellous  bone  as  a whole  into  a 
problem  of  determining  the  remodeling  properties  of  a single  trabecular  plate. 

Unfortunately,  since  bone  remodeling  can  only  occur  with  the  action  of 
living  bone  cells,  the  measurement  of  theRemodeling  Function  can  only  be 
accomplished  with  in-vivo  testing.  At  the  present  time,  the  most  likely 
candidate  as  a site  for  in-vivo  testing  is  the  calcaneus  (or  heel  bone ) of  a 
primate.  This  experimental  site  is  desirable  for  the  following  reasons: 

(1)  the  calcaneus  has  distinctive  trabecular  patterns  which  show 
up  very  clearly  on  radiographs,  and 

(2)  the  loading  patterns  and  geometric  configuration  of  a 
calcaneus  are  such  that  the  adaptive  elastic  problem  may  be 
considered  symmetric  in  the  medial-anterior  plane. 

Figure  18  illustrates  an  experimental  test  configuration  for  measuring  the 
Remodeling  Function  of  the  calcaneus.  One  foot  is  immobilized  from  the  ankle 
down  and  a calibrated  load  is  applied  to  the  heel,  while  the  other  foot  is  used  as 
a control.  Periodically  radiographs  can  be  taken  and  an  autoniciic  polarimeter 
used  to  determine  the  strain  in  each  region  of  the  calcaneous  as  a function  of 
the  calibrated  load.  Observing  the  rate  change  in  the  trabecular  adaptation 
tensor  M as  a function  of  the  strain  and  the  A 
adaptation,  one  can  then  evaluate  the  function  R..(E,M). 
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IV.  4 CLOSURE 


The  Theory  of  Adaptive  Elasticity  presented  here  should  prove  to  be  a 
useful  tool  in  the  practice  of  orthopedics.  The  governing  equations  of 
Adaptive  Elasticity  not  only  model  the  important  physiological  phenomena 
of  stres 9 -adaptation  of  bone,  but  they  also  might  enable  a clearer  understand- 
ing of  the  etiology  of  certain  muscular-skeletal  disorders  such  as  osteoporosis 
and  scoliosis.  As  an  added  bonus,  our  mathematical  analysis  has  suggested 
the  design  oi  new  instrumentation,  an  automatic  polarimeter,  which  can 
quantitatively  measure  the  trabecular  adaptation  of  cancellous  bone  from  a 
radiograph.  This  instrumentation  would  have  immediate  clinical  applications 
in  the  diagnosis  and  treatment  of  osteoporosis. 
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